<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">

<head>

<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />

<meta name="viewport" content="width=device-width, initial-scale=1">

<meta name="author" content="Peter Bourke, Geert van Geest" />

<meta name="date" content="2018-02-13" />

<title>polymapR: linkage mapping in outcrossing polyploids</title>



<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
div.sourceCode { overflow-x: auto; }
table.sourceCode, tr.sourceCode, td.lineNumbers, td.sourceCode {
  margin: 0; padding: 0; vertical-align: baseline; border: none; }
table.sourceCode { width: 100%; line-height: 100%; }
td.lineNumbers { text-align: right; padding-right: 4px; padding-left: 4px; color: #aaaaaa; border-right: 1px solid #aaaaaa; }
td.sourceCode { padding-left: 5px; }
code > span.kw { color: #007020; font-weight: bold; } /* Keyword */
code > span.dt { color: #902000; } /* DataType */
code > span.dv { color: #40a070; } /* DecVal */
code > span.bn { color: #40a070; } /* BaseN */
code > span.fl { color: #40a070; } /* Float */
code > span.ch { color: #4070a0; } /* Char */
code > span.st { color: #4070a0; } /* String */
code > span.co { color: #60a0b0; font-style: italic; } /* Comment */
code > span.ot { color: #007020; } /* Other */
code > span.al { color: #ff0000; font-weight: bold; } /* Alert */
code > span.fu { color: #06287e; } /* Function */
code > span.er { color: #ff0000; font-weight: bold; } /* Error */
code > span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
code > span.cn { color: #880000; } /* Constant */
code > span.sc { color: #4070a0; } /* SpecialChar */
code > span.vs { color: #4070a0; } /* VerbatimString */
code > span.ss { color: #bb6688; } /* SpecialString */
code > span.im { } /* Import */
code > span.va { color: #19177c; } /* Variable */
code > span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code > span.op { color: #666666; } /* Operator */
code > span.bu { } /* BuiltIn */
code > span.ex { } /* Extension */
code > span.pp { color: #bc7a00; } /* Preprocessor */
code > span.at { color: #7d9029; } /* Attribute */
code > span.do { color: #ba2121; font-style: italic; } /* Documentation */
code > span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code > span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code > span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
</style>



<link href="data:text/css;charset=utf-8,body%20%7B%0Abackground%2Dcolor%3A%20%23fff%3B%0Amargin%3A%201em%20auto%3B%0Amax%2Dwidth%3A%20700px%3B%0Aoverflow%3A%20visible%3B%0Apadding%2Dleft%3A%202em%3B%0Apadding%2Dright%3A%202em%3B%0Afont%2Dfamily%3A%20%22Open%20Sans%22%2C%20%22Helvetica%20Neue%22%2C%20Helvetica%2C%20Arial%2C%20sans%2Dserif%3B%0Afont%2Dsize%3A%2014px%3B%0Aline%2Dheight%3A%201%2E35%3B%0A%7D%0A%23header%20%7B%0Atext%2Dalign%3A%20center%3B%0A%7D%0A%23TOC%20%7B%0Aclear%3A%20both%3B%0Amargin%3A%200%200%2010px%2010px%3B%0Apadding%3A%204px%3B%0Awidth%3A%20400px%3B%0Aborder%3A%201px%20solid%20%23CCCCCC%3B%0Aborder%2Dradius%3A%205px%3B%0Abackground%2Dcolor%3A%20%23f6f6f6%3B%0Afont%2Dsize%3A%2013px%3B%0Aline%2Dheight%3A%201%2E3%3B%0A%7D%0A%23TOC%20%2Etoctitle%20%7B%0Afont%2Dweight%3A%20bold%3B%0Afont%2Dsize%3A%2015px%3B%0Amargin%2Dleft%3A%205px%3B%0A%7D%0A%23TOC%20ul%20%7B%0Apadding%2Dleft%3A%2040px%3B%0Amargin%2Dleft%3A%20%2D1%2E5em%3B%0Amargin%2Dtop%3A%205px%3B%0Amargin%2Dbottom%3A%205px%3B%0A%7D%0A%23TOC%20ul%20ul%20%7B%0Amargin%2Dleft%3A%20%2D2em%3B%0A%7D%0A%23TOC%20li%20%7B%0Aline%2Dheight%3A%2016px%3B%0A%7D%0Atable%20%7B%0Amargin%3A%201em%20auto%3B%0Aborder%2Dwidth%3A%201px%3B%0Aborder%2Dcolor%3A%20%23DDDDDD%3B%0Aborder%2Dstyle%3A%20outset%3B%0Aborder%2Dcollapse%3A%20collapse%3B%0A%7D%0Atable%20th%20%7B%0Aborder%2Dwidth%3A%202px%3B%0Apadding%3A%205px%3B%0Aborder%2Dstyle%3A%20inset%3B%0A%7D%0Atable%20td%20%7B%0Aborder%2Dwidth%3A%201px%3B%0Aborder%2Dstyle%3A%20inset%3B%0Aline%2Dheight%3A%2018px%3B%0Apadding%3A%205px%205px%3B%0A%7D%0Atable%2C%20table%20th%2C%20table%20td%20%7B%0Aborder%2Dleft%2Dstyle%3A%20none%3B%0Aborder%2Dright%2Dstyle%3A%20none%3B%0A%7D%0Atable%20thead%2C%20table%20tr%2Eeven%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0A%7D%0Ap%20%7B%0Amargin%3A%200%2E5em%200%3B%0A%7D%0Ablockquote%20%7B%0Abackground%2Dcolor%3A%20%23f6f6f6%3B%0Apadding%3A%200%2E25em%200%2E75em%3B%0A%7D%0Ahr%20%7B%0Aborder%2Dstyle%3A%20solid%3B%0Aborder%3A%20none%3B%0Aborder%2Dtop%3A%201px%20solid%20%23777%3B%0Amargin%3A%2028px%200%3B%0A%7D%0Adl%20%7B%0Amargin%2Dleft%3A%200%3B%0A%7D%0Adl%20dd%20%7B%0Amargin%2Dbottom%3A%2013px%3B%0Amargin%2Dleft%3A%2013px%3B%0A%7D%0Adl%20dt%20%7B%0Afont%2Dweight%3A%20bold%3B%0A%7D%0Aul%20%7B%0Amargin%2Dtop%3A%200%3B%0A%7D%0Aul%20li%20%7B%0Alist%2Dstyle%3A%20circle%20outside%3B%0A%7D%0Aul%20ul%20%7B%0Amargin%2Dbottom%3A%200%3B%0A%7D%0Apre%2C%20code%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0Aborder%2Dradius%3A%203px%3B%0Acolor%3A%20%23333%3B%0Awhite%2Dspace%3A%20pre%2Dwrap%3B%20%0A%7D%0Apre%20%7B%0Aborder%2Dradius%3A%203px%3B%0Amargin%3A%205px%200px%2010px%200px%3B%0Apadding%3A%2010px%3B%0A%7D%0Apre%3Anot%28%5Bclass%5D%29%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0A%7D%0Acode%20%7B%0Afont%2Dfamily%3A%20Consolas%2C%20Monaco%2C%20%27Courier%20New%27%2C%20monospace%3B%0Afont%2Dsize%3A%2085%25%3B%0A%7D%0Ap%20%3E%20code%2C%20li%20%3E%20code%20%7B%0Apadding%3A%202px%200px%3B%0A%7D%0Adiv%2Efigure%20%7B%0Atext%2Dalign%3A%20center%3B%0A%7D%0Aimg%20%7B%0Abackground%2Dcolor%3A%20%23FFFFFF%3B%0Apadding%3A%202px%3B%0Aborder%3A%201px%20solid%20%23DDDDDD%3B%0Aborder%2Dradius%3A%203px%3B%0Aborder%3A%201px%20solid%20%23CCCCCC%3B%0Amargin%3A%200%205px%3B%0A%7D%0Ah1%20%7B%0Amargin%2Dtop%3A%200%3B%0Afont%2Dsize%3A%2035px%3B%0Aline%2Dheight%3A%2040px%3B%0A%7D%0Ah2%20%7B%0Aborder%2Dbottom%3A%204px%20solid%20%23f7f7f7%3B%0Apadding%2Dtop%3A%2010px%3B%0Apadding%2Dbottom%3A%202px%3B%0Afont%2Dsize%3A%20145%25%3B%0A%7D%0Ah3%20%7B%0Aborder%2Dbottom%3A%202px%20solid%20%23f7f7f7%3B%0Apadding%2Dtop%3A%2010px%3B%0Afont%2Dsize%3A%20120%25%3B%0A%7D%0Ah4%20%7B%0Aborder%2Dbottom%3A%201px%20solid%20%23f7f7f7%3B%0Amargin%2Dleft%3A%208px%3B%0Afont%2Dsize%3A%20105%25%3B%0A%7D%0Ah5%2C%20h6%20%7B%0Aborder%2Dbottom%3A%201px%20solid%20%23ccc%3B%0Afont%2Dsize%3A%20105%25%3B%0A%7D%0Aa%20%7B%0Acolor%3A%20%230033dd%3B%0Atext%2Ddecoration%3A%20none%3B%0A%7D%0Aa%3Ahover%20%7B%0Acolor%3A%20%236666ff%3B%20%7D%0Aa%3Avisited%20%7B%0Acolor%3A%20%23800080%3B%20%7D%0Aa%3Avisited%3Ahover%20%7B%0Acolor%3A%20%23BB00BB%3B%20%7D%0Aa%5Bhref%5E%3D%22http%3A%22%5D%20%7B%0Atext%2Ddecoration%3A%20underline%3B%20%7D%0Aa%5Bhref%5E%3D%22https%3A%22%5D%20%7B%0Atext%2Ddecoration%3A%20underline%3B%20%7D%0A%0Acode%20%3E%20span%2Ekw%20%7B%20color%3A%20%23555%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%0Acode%20%3E%20span%2Edt%20%7B%20color%3A%20%23902000%3B%20%7D%20%0Acode%20%3E%20span%2Edv%20%7B%20color%3A%20%2340a070%3B%20%7D%20%0Acode%20%3E%20span%2Ebn%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Efl%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Ech%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Est%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Eco%20%7B%20color%3A%20%23888888%3B%20font%2Dstyle%3A%20italic%3B%20%7D%20%0Acode%20%3E%20span%2Eot%20%7B%20color%3A%20%23007020%3B%20%7D%20%0Acode%20%3E%20span%2Eal%20%7B%20color%3A%20%23ff0000%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%0Acode%20%3E%20span%2Efu%20%7B%20color%3A%20%23900%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%20code%20%3E%20span%2Eer%20%7B%20color%3A%20%23a61717%3B%20background%2Dcolor%3A%20%23e3d2d2%3B%20%7D%20%0A" rel="stylesheet" type="text/css" />

</head>

<body>




<h1 class="title toc-ignore">polymapR: linkage mapping in outcrossing polyploids</h1>
<h4 class="author"><em>Peter Bourke, Geert van Geest</em></h4>
<h4 class="date"><em>2018-02-13</em></h4>



<p>In this tutorial, we will go step by step through the basic steps of performing linkage analysis in polyploids using <code>polymapR</code>. Most of the theory behind the functions described here can be found in a bioRXiv preprint manuscript of Bourke et al. (2017) at <a href="https://doi.org/10.1101/228817" class="uri">https://doi.org/10.1101/228817</a>. The example dataset is a simulated dataset of an hypothetical tetraploid outcrossing plant species with five chromosomes. Chromosomes have random pairing behaviour. The population consists of 207 individuals originating from an F1 cross between two non-inbred parents. Frequently occurring data imperfections like missing values, genotyping errors and duplicated individuals are taken into account.</p>
<div id="table-of-contents" class="section level2">
<h2>Table of Contents</h2>
<ol style="list-style-type: decimal">
<li><a href="#sect1">A (very) short introduction to R and the <code>polymapR</code> package.</a></li>
<li><a href="#sect2">Install <code>polymapR</code></a></li>
<li><a href="#sect3">Logging function calls and output</a></li>
<li><a href="#sect4">Data importing</a></li>
<li><a href="#sect5">Generate summary data</a></li>
<li><a href="#sect6">Convert marker dosages to simple segregations; remove non-segregating data</a></li>
<li><a href="#sect7">Quality checks on marker data</a></li>
<li><a href="#sect8">Simplex x nulliplex markers – defining chromosomes and homologues</a></li>
<li><a href="#sect9">Assigning SxS and DxN markers and consensus linkage group (LG) names</a></li>
<li><a href="#sect10">Assign all other markertypes</a></li>
<li><a href="#sect11">Finish the linkage analysis</a></li>
<li><a href="#sect12"><em>Optional</em>: Marker binning</a></li>
<li><a href="#sect13"><em>Optional</em>: Writing to a .pwd file</a></li>
<li><a href="#sect14">Marker ordering</a></li>
<li><a href="#sect15">Plotting a map</a></li>
<li><a href="#sect16">Evaluating map quality</a></li>
<li><a href="#sect17">Preferential pairing</a></li>
<li><a href="#sect18">Concluding remarks</a></li>
<li><a href="#sect19">References</a></li>
</ol>
<p><a name="sect1"></a></p>
</div>
<div id="a-very-short-introduction-to-r-and-the-polymapr-package." class="section level2">
<h2>1. A (very) short introduction to R and the <code>polymapR</code> package.</h2>
<p>In R, a function can be called by typing <code>function_name()</code> in the console and pressing <code>[Enter]</code>. (Optional) arguments go within the parentheses, followed by a <code>=</code> sign and the value you want to give to the argument. Arguments are separated by a comma. All possible arguments and default values of a function are in the help file, which is called by typing <code>?</code> followed by the function name. Getting the help file for the commonly used function <code>seq</code> looks like this:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">?seq</code></pre></div>
<p>And generating a sequence of numbers from 2 to 14 by an increment of 3 using <code>seq</code> looks like this:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">seq</span>(<span class="dt">from =</span> <span class="dv">2</span>,
    <span class="dt">to =</span> <span class="dv">14</span>,
    <span class="dt">by =</span> <span class="dv">3</span>)</code></pre></div>
<pre><code>## [1]  2  5  8 11 14</code></pre>
<p>Although copy pasting the commands from this document to your R console will do for this tutorial, we recommend to have a bit more experience in using R than above before using <code>polymapR</code>. This will allow you better to troubleshoot and check out data output. A good place to start would be <a href="http://tryr.codeschool.com" class="uri">http://tryr.codeschool.com</a> or any of the plethora of R tutorials on the web.</p>
<p>The <code>polymapR</code> package is literally a collection of functions, each defined by their function name and arguments. It relies on dosage data of molecular markers, and therefore is logically used after the packages <code>fitTetra</code> or <code>fitPoly</code>. As for any project, it is best to first create a new folder which will be the main working directory for the mapping project. Set the working directory to your new folder (<code>setwd()</code> or use the drop-down menu) and copy the marker dosage file from <code>fitTetra</code> there.</p>
<p><a name="sect2"></a></p>
</div>
<div id="install-polymapr" class="section level2">
<h2>2. Install <code>polymapR</code></h2>
<p><code>polymapR</code> is available on CRAN at <a href="https://cran.R-project.org/package=polymapR" class="uri">https://cran.R-project.org/package=polymapR</a> which is probably where you accessed this vignette. As with other packages on CRAN, it can be installed directly from within R using the command:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">install.packages</span>(<span class="st">&quot;polymapR&quot;</span>)</code></pre></div>
<p><a name="sect3"></a></p>
</div>
<div id="logging-function-calls-and-output" class="section level2">
<h2>3. Logging function calls and output</h2>
<p>Most exported function in <code>polymapR</code> have a <code>log</code> argument. The default setting is <code>NULL</code>. In that case, all messages are returned to the standard output, which means your R console if you are calling the functions interactively. However, if it is set to a file name (e.g. <code>log.Rmd</code> or <code>log.txt</code>) all function calls and output are written to that file. It automatically creates one if it does not exist yet. The log file is written in R markdown (<a href="http://rmarkdown.rstudio.com/" class="uri">http://rmarkdown.rstudio.com/</a>). The advantage of markdown is that it is readable as plain text, but it can also easily be rendered into a nicely formatted html or pdf. In order to do so, give your logfile a .Rmd extension and use the <code>knitr</code> package to turn it into a pdf, docx or html.</p>
<p><a name="sect4"></a></p>
</div>
<div id="data-importing" class="section level2">
<h2>4. Data importing</h2>
<div id="reading-in-marker-dosage-data" class="section level3">
<h3>4.1 Reading in marker dosage data</h3>
<p>We assume that the user has used <code>fitTetra</code> or some alternative method to convert signal intensity ratios into marker dosages. <code>polymapR</code> takes as its starting point such data, but we will impose some further quality checks before we can begin mapping (such as allowed numbers of missing values). We also need to check the markers for skewness (differences in the observed and expected marker segregation) which can cause problems later in the clustering. For this, we have incorporated a number of functions developed for the fitTetra package <span class="citation">[1]</span>, and included in the currently-unpublished R package <code>fitPolyTools</code> developed by Roeland Voorrips.</p>
<blockquote>
<h6 id="but-first-a-note-on-polyploid-terminology">But first, a note on polyploid terminology</h6>
<p>In the <code>polymapR</code> package, and throughout this vignette, we often refer to terms like “simplex”, “duplex” or “nulliplex” etc. These terms describe the dosage scores of SNP markers, essentially the count of the number of “alternative” alleles present at a particular locus in an individual. A simplex x nulliplex marker is therefore a marker that for which the maternal dosage is 1 (mothers always come first) and paternal dosage is 0. To save space, we also refer to these markers as SxN markers, which is synonymous with 1x0 markers. Sometimes statements are made concerning simplex x nulliplex markers - it should be understood that these statements also usually apply to nulliplex x simplex markers - markers for which the segregating allele is carried by the father rather than the mother.</p>
</blockquote>
<p>The first step is to import the marker data. The layout of the dataset should be marker name, parental scores followed by offspring scores (so-called wide data format):</p>
<table>
<thead>
<tr class="header">
<th align="left">Marker</th>
<th align="right">P1</th>
<th align="right">P2</th>
<th align="right">F1.001</th>
<th align="right">F1.002..</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td align="left">mymarker1</td>
<td align="right">1</td>
<td align="right">0</td>
<td align="right">1</td>
<td align="right">1</td>
</tr>
<tr class="even">
<td align="left">mymarker2</td>
<td align="right">2</td>
<td align="right">4</td>
<td align="right">3</td>
<td align="right">2</td>
</tr>
<tr class="odd">
<td align="left">mymarker3</td>
<td align="right">0</td>
<td align="right">3</td>
<td align="right">2</td>
<td align="right">2</td>
</tr>
</tbody>
</table>
<p>In cases where the parents have been genotyped more than once, a consensus parental score is required. Missing parental scores are not allowed at this stage, and should have been imputed already if desired.</p>
<p>The data should be in some format that can be read easily into R, such as .csv, .dat or .txt file. In the example a .csv file is read in. For other file types check out <code>?read.delim</code>.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ALL_dosages &lt;-<span class="st"> </span><span class="kw">read.csv</span>(<span class="st">&quot;tetraploid_dosages.csv&quot;</span>,
                        <span class="dt">stringsAsFactors =</span> <span class="ot">FALSE</span>,
                        <span class="dt">row.names =</span> <span class="dv">1</span>) <span class="co">#first column contains rownames</span></code></pre></div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">class</span>(ALL_dosages)</code></pre></div>
<pre><code>## [1] &quot;data.frame&quot;</code></pre>
<p>By default, <code>read.csv</code> returns data in a data frame. All <code>polymapR</code> functions that use dosage data only accept dosage data when the dosages are in a matrix. A matrix is a two-dimensional R object with all elements of the same type (e.g. character, integer, numeric). Markernames and individual names should be specified in rownames and columnnames respectively. To use the data with functions, the data should be converted from a data.frame in to a matrix by the <code>as.matrix</code> function:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ALL_dosages &lt;-<span class="st"> </span><span class="kw">as.matrix</span>(ALL_dosages)

<span class="kw">class</span>(ALL_dosages)</code></pre></div>
<pre><code>## [1] &quot;matrix&quot;</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">head</span>(ALL_dosages[,<span class="dv">1</span>:<span class="dv">5</span>])</code></pre></div>
<pre><code>##             P1 P2 F1_001 F1_002 F1_003
## Zm_rs005505  2  2      2      3      2
## Zm_ts097822  1  2      1      1      2
## Ac_ws072452  1  0      1     NA      1
## Ac_ts073123  0  2      1     NA      2
## Ap_ws071152  3  1      3      2     NA
## St_rs076767  1  0      1      1      1</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">dim</span>(ALL_dosages)</code></pre></div>
<pre><code>## [1] 3000  209</code></pre>
<p>If you are testing the package and would like a sample dataset to work with, a set of tetraploid marker dosage data is provided within the package itself. First load the package and then the dosage data as follows:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(polymapR)
<span class="kw">data</span>(ALL_dosages)</code></pre></div>
</div>
<div id="checking-for-skewness" class="section level3">
<h3>4.2 Checking for skewness</h3>
<p>The next step before we proceed further is to check whether the marker scores in the F1 correspond to those expected according to the parental dosages. We do this using the <code>checkF1</code> function (from <code>fitPolyTools</code>). However, we first load the <code>polymapR</code> package into the global environment:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(polymapR)</code></pre></div>
<p>Now run <code>checkF1</code>:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">F1checked &lt;-<span class="st"> </span><span class="kw">checkF1</span>(<span class="dt">dosage_matrix =</span> ALL_dosages,<span class="dt">parent1 =</span> <span class="st">&quot;P1&quot;</span>,<span class="dt">parent2 =</span> <span class="st">&quot;P2&quot;</span>,
                     <span class="dt">F1 =</span> <span class="kw">colnames</span>(ALL_dosages)[<span class="dv">3</span>:<span class="kw">ncol</span>(ALL_dosages)],
                     <span class="dt">polysomic =</span> <span class="ot">TRUE</span>, <span class="dt">disomic =</span> <span class="ot">FALSE</span>, <span class="dt">mixed =</span> <span class="ot">FALSE</span>, <span class="dt">ploidy =</span> <span class="dv">4</span>)

<span class="kw">head</span>(F1checked)</code></pre></div>
<pre><code>##   m  MarkerName parent1 parent2 F1_0 F1_1 F1_2 F1_3 F1_4 F1_NA P1 P2
## 1 1 Zm_rs005505       2       2    6   39   90   45    4    23  2  2
## 2 2 Zm_ts097822       1       2   14   79   83   10    0    21  1  2
## 3 3 Ac_ws072452       1       0  102   83    0    0    0    22  1  0
## 4 4 Ac_ts073123       0       2   27  121   37    0    0    22  0  2
## 5 5 Ap_ws071152       3       1    1   48   77   55    0    26  3  1
## 6 6 St_rs076767       1       0   92   88    0    0    0    27  1  0
##   bestfit frqInvalid_bestfit Pvalue_bestfit matchParent_bestfit
## 1 18I81_0             0.0000         0.9187                 Yes
## 2  1551_0             0.0000         0.4724                 Yes
## 3    11_0             0.0000         0.1624                 Yes
## 4   141_0             0.0000         0.4160                 Yes
## 5   121_1             0.0055         0.1165                 Yes
## 6    11_0             0.0000         0.7656                 Yes
##   bestParentfit frqInvalid_bestParentfit Pvalue_bestParentfit
## 1       18I81_0                   0.0000               0.9187
## 2        1551_0                   0.0000               0.4724
## 3          11_0                   0.0000               0.1624
## 4         141_0                   0.0000               0.4160
## 5         121_1                   0.0055               0.1165
## 6          11_0                   0.0000               0.7656
##   matchParent_bestParentfit q1_segtypefit q2_parents q3_fracscored
## 1                       Yes        1.0000          1        0.7222
## 2                       Yes        1.0000          1        0.7464
## 3                       Yes        1.0000          1        0.7343
## 4                       Yes        1.0000          1        0.7343
## 5                       Yes        0.9278          1        0.6860
## 6                       Yes        1.0000          1        0.6739
##   qall_mult qall_weights
## 1    0.7222       0.9383
## 2    0.7464       0.9436
## 3    0.7343       0.9410
## 4    0.7343       0.9410
## 5    0.6365       0.8901
## 6    0.6739       0.9275</code></pre>
<p>There is the possibility to check for shifted markers (incorrectly genotyped) but we will not concern ourselves with that right now. The main arguments to be specified are what the parental identifiers are (<code>parent1</code> and <code>parent2</code>), what the identifiers of the F1 population are (argument <code>F1</code> - note we just take the column names of the dosages matrix from column 3 to the end), and which sort of inheritance model we wish to compare against (here we have specified <code>polysomic</code> to be <code>TRUE</code> and both <code>disomic</code> and <code>mixed</code> to be <code>FALSE</code>. Note that <code>mixed</code> allows for one parent to be fully disomic and the other fully polysomic, but does not refer to the situation of a mixed inheritance pattern within one parent, associated with segmental allopolyploidy - this is because we have no fixed expectation to compare with if this is the case). The argument <code>ploidy</code> specifies the ploidy level of parent 1, if parent 2 has a different ploidy (for instance in a tetraploid x diploid cross) then we can use the <code>ploidy2</code> argument. The results are stored in the data.frame <code>F1checked</code>, which has quite a few columns. We are interested in seeing whether there is consistency between the parental scores and the F1. In general, if there is good consistency between the parental and offspring scores we will have a high value for <code>qall_mult</code> and <code>qall_weights</code>, which are aggregated quality measures of <code>q1</code>, <code>q2</code> and <code>q3</code>. In the sample datatset provided with <code>polymapR</code> there are no skewed markers present. However, you may find it useful to examine the output of <code>checkF1</code> and remove markers which show conflict between the parental and offspring scores, as these markers will likely cause issues in the later mapping steps.</p>
<p><a name="sect5"></a></p>
</div>
</div>
<div id="generate-summary-data" class="section level2">
<h2>5. Generate summary data</h2>
<p>Now we summarise the marker data in terms of the different segregation types, number of missing values, number of inadmissible scores (which could be genotyping errors or possible double reduction scores). Run the summary function on the marker data as follows:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mds &lt;-<span class="st"> </span><span class="kw">marker_data_summary</span>(<span class="dt">dosage_matrix =</span> ALL_dosages,
                           <span class="dt">ploidy =</span> <span class="dv">4</span>,
                           <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>,
                           <span class="dt">parent1 =</span> <span class="st">&quot;P1&quot;</span>,
                           <span class="dt">parent2 =</span> <span class="st">&quot;P2&quot;</span>,
                           <span class="dt">progeny_incompat_cutoff =</span> <span class="fl">0.05</span>)</code></pre></div>
<pre><code>## Calculating parental info...
## Checking compatability between parental and offspring scores...
## 
## ####parental_info
## 
##         P2_0   P2_1   P2_2   P2_3   P2_4
## -----  -----  -----  -----  -----  -----
## P1_0       0    312    227     49      3
## P1_1     300    625    356     67      4
## P1_2     230    323    175     29      2
## P1_3      66     65     30      5      0
## P1_4       2      5      3      0    122
## 
## ####offspring_incompatible
## 
##         P2_0   P2_1   P2_2   P2_3   P2_4
## -----  -----  -----  -----  -----  -----
## P1_0      NA   0.00   0.00   0.56   2.42
## P1_1    0.00   0.00   0.00   0.28   0.85
## P1_2    0.00   0.00   0.00   0.12   0.00
## P1_3    0.45   0.33   0.14   0.19     NA
## P1_4    0.97   0.48   0.64     NA   1.07
## 
## ####Incompatible individuals:
## 
## None</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">pq_before_convert &lt;-<span class="st"> </span><span class="kw">parental_quantities</span>(<span class="dt">dosage_matrix =</span> ALL_dosages, 
                                         <span class="dt">las =</span> <span class="dv">2</span>)</code></pre></div>
<p><img src="" /><!-- --></p>
<pre><code>##              1x0   2x0   3x0   4x0   0x1   1x1   2x1   3x1   4x1   0x2   1x2   2x2   3x2   4x2   0x3   1x3   2x3   3x3   0x4   1x4   2x4   4x4
## ----------  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----  ----
## frequency    300   230    66     2   312   625   323    65     5   227   356   175    30     3    49    67    29     5     3     4     2   122</code></pre>
<p>If you want to know exactly what these functions are doing and what can be typed as an argument just use <code>?marker_data_summary</code> and/or <code>?parental_quantities</code>.</p>
<p><a name="sect6"></a></p>
</div>
<div id="convert-marker-dosages-to-simple-segregations-remove-non-segregating-data" class="section level2">
<h2>6. Convert marker dosages to simple segregations; remove non-segregating data</h2>
<p>Simplex x nulliplex markers segregate in a 1:1 ratio in the offspring, but this segregation is also possible with triplex x nulliplex, triplex x quadriplex and simplex x quadriplex, and the segregating allele in each case comes from the same parent. We can therefore re-code all of these segregation types into the simplest 1 x 0 case, taking care that we handle possible double reduction scores correctly. A list of the conversions in tetraploid we apply which simplify the subsequent analysis and mapping work are:</p>
<table>
<tbody>
<tr class="odd">
<td align="left">3x0,3x4,1x4 -&gt; 1x0</td>
</tr>
<tr class="even">
<td align="left">0x3,4x3,4x1 -&gt; 0x1</td>
</tr>
<tr class="odd">
<td align="left">2x4 -&gt; 2x0</td>
</tr>
<tr class="even">
<td align="left">4x2 -&gt; 0x2</td>
</tr>
<tr class="odd">
<td align="left">3x3 -&gt; 1x1</td>
</tr>
<tr class="even">
<td align="left">3x1 -&gt; 1x3</td>
</tr>
<tr class="odd">
<td align="left">3x2 -&gt; 1x2</td>
</tr>
<tr class="even">
<td align="left">2x3 -&gt; 2x1</td>
</tr>
</tbody>
</table>
<p>Note that some classes such as 2x2 cannot be converted into any other category, and we do not convert non-segregating classes such as 4x4 as these are of no further use for mapping or QTL analysis. To convert the marker data:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">segregating_data &lt;-<span class="st"> </span><span class="kw">convert_marker_dosages</span>(<span class="dt">dosage_matrix =</span> ALL_dosages)</code></pre></div>
<pre><code>## Warning in convert_marker_dosages(dosage_matrix = ALL_dosages): There are dosages greater than 4 or less than 0 in the dataset.
##             If they include parental dosages, markers are removed from the output.
##             Otherwise the dosage is made missing.</code></pre>
<pre><code>## 
## ####Marker dosage frequencies:
## 
##         P2_0   P2_1   P2_2   P2_3
## -----  -----  -----  -----  -----
## P1_0       0    366    230      0
## P1_1     370    630    386    132
## P1_2     232    352    175      0
## 
## markers not converted: 2615
## 
## markers 1 parent converted: 193
## 
## markers 2 parents converted: 192
## 
## non-segregating markers deleted: 127</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">pq_after_convert &lt;-<span class="st"> </span><span class="kw">parental_quantities</span>(<span class="dt">dosage_matrix =</span> segregating_data)</code></pre></div>
<p><img src="" /><!-- --></p>
<pre><code>##              1x0   2x0   0x1   1x1   2x1   0x2   1x2   2x2   1x3
## ----------  ----  ----  ----  ----  ----  ----  ----  ----  ----
## frequency    370   232   366   630   352   230   386   175   132</code></pre>
<p><a name="sect7"></a></p>
</div>
<div id="quality-checks-on-marker-data" class="section level2">
<h2>7. Quality checks on marker data</h2>
<p>There are a number of quality checks which are needed before we begin mapping. For example, we do not want to include markers which have too many missing values (“NA” scores) as this may be an indication that the marker has performed poorly and may have many errors associated with it. We would also like to screen out individuals that have many missing values as this might be an indication that the DNA quality of that sample was poor, or duplicate individuals which give no extra information and should probably also be removed / combined. Because there are no fixed rules for how stringent we should be with quality checks (5% missing values ok? 10% missing values ok? 20%?), some visual aids are produced in order to help with these decisions.</p>
<div id="missing-value-rate-per-marker" class="section level3">
<h3>7.1 Missing value rate per marker</h3>
<p>We first screen the marker data for missing values using <code>screen_for_NA_values</code>. This function allows the user to remove rows or columns with a certain rate of missing values from the dataset. The dosage matrix, which is specified as <code>segregating_data</code> contains markers in rows and individuals in columns. In R, rows are usually specified by the number 1, and columns by 2. The argument <code>margin</code> lets you specify by margin number whether you want to screen markers (<code>margin = 1</code>) or individuals (<code>margin = 2</code>) If the argument <code>cutoff</code> is set to <code>NULL</code>, it takes user input on the rate of markers to be screened. In the example below <code>cutoff</code> is specified, but because decisions are made after data inspection, we recommend to decide on the cutoff value after inspection of the data by specifying <code>cutoff = NULL</code>. A threshold of 10% missing values is sensible, as a compromise between high-confidence marker data versus keeping a large proportion of the markers for mapping and QTL analysis. In general, the most problematic markers tend to have higher rates of missing values so any reasonable threshold here should be fine. If problems are encountered later in the mapping, it might be worthwhile to re-run this step with a more stringent (e.g. 5%) threshold.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">screened_data &lt;-<span class="st"> </span><span class="kw">screen_for_NA_values</span>(<span class="dt">dosage_matrix =</span> segregating_data, 
                                      <span class="dt">margin =</span> <span class="dv">1</span>, <span class="co"># margin 1 means markers</span>
                                      <span class="dt">cutoff =</span>  <span class="fl">0.10</span>,
                                      <span class="dt">print.removed =</span> <span class="ot">FALSE</span>) </code></pre></div>
<p><img src="" /><!-- --></p>
<pre><code>## 
## Number of removed markers :  1458 
## 
## 
## There are now 1415 markers leftover in the dataset.</code></pre>
</div>
<div id="missing-value-rate-per-individual" class="section level3">
<h3>7.2 Missing value rate per individual</h3>
<p>As in step 7.1 we can also visualise the rate of missing values per individual, which can help us to make a decision about whether certain individuals should be removed from the dataset before commencing the mapping. In order to screen the individuals for missing values we enter:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">screened_data2 &lt;-<span class="st"> </span><span class="kw">screen_for_NA_values</span>(<span class="dt">dosage_matrix =</span> screened_data, 
                                       <span class="dt">cutoff =</span> <span class="fl">0.1</span>, 
                                       <span class="dt">margin =</span> <span class="dv">2</span>, <span class="co">#margin = 2 means columns</span>
                                       <span class="dt">print.removed =</span> <span class="ot">FALSE</span>)</code></pre></div>
<p><img src="" /><!-- --></p>
<pre><code>## 
## Number of removed individuals :  2 
## 
## 
## There are now 205 individuals leftover in the dataset.</code></pre>
<p>Again, it is not clear what an acceptable threshold for the rate of missing values should be, but we recommend a cut-off rate of 0.1 – 0.15 (10-15%) missing values allowed.</p>
</div>
<div id="duplicate-individuals" class="section level3">
<h3>7.3 Duplicate individuals</h3>
<p>It is also desirable to be able to check whether there are any duplicated individuals in the population as can sometimes happen (or due to mix-ups in the DNA preparation and genotyping). The function <code>screen_for_duplicate_individuals</code> does this for you. Note that a duplicate can remain in the dataset, but they add no further information for mapping or QTL analysis, and may even distort the analysis (by giving an inaccurate population size). As for the function <code>screen_for_NA_values</code>, cutoff can be set to <code>NULL</code> to accept user input. Especially for <code>screen_for_duplicate_individuals</code> it makes sense to make a decision on the cut off after inspection of the data and figures. To screen the data for duplicate individuals, enter:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">screened_data3 &lt;-<span class="st"> </span><span class="kw">screen_for_duplicate_individuals</span>(<span class="dt">dosage_matrix =</span> screened_data2, 
                                                   <span class="dt">cutoff =</span> <span class="fl">0.95</span>, 
                                                   <span class="dt">plot_cor =</span> T)</code></pre></div>
<p><img src="" /><!-- --></p>
<pre><code>## 
## Combining F1_060 &amp; F1_081 into F1_060
## 
## Combining F1_046 &amp; F1_085 into F1_046
## 
## Combining F1_032 &amp; F1_100 into F1_032
## 
## Combining F1_079 &amp; F1_101 into F1_079
## 
## Combining F1_106 &amp; F1_113 into F1_106
## 
## Combining F1_112 &amp; F1_199 into F1_112
## 
## #### 6 individuals removed:
## 
##      _        _        _        _       
## [1,] &quot;F1_081&quot; &quot;F1_085&quot; &quot;F1_100&quot; &quot;F1_114&quot;
## [2,] &quot;F1_101&quot; &quot;F1_113&quot; &quot;F1_199&quot; &quot;&quot;      
## _        _        _        _      
## -------  -------  -------  -------
## F1_081   F1_085   F1_100   F1_114 
## F1_101   F1_113   F1_199</code></pre>
<p>Here we see that there were a number of duplicate pairs of individuals identified (red points in the output graph). Each duplicate pair is merged into a consensus individual (conflicts are made missing), keeping the name of the first individual. This is summarised in the output as well.</p>
</div>
<div id="duplicated-markers" class="section level3">
<h3>7.4 Duplicated markers</h3>
<p>Finally, we can also check for duplicated markers (those which give us no extra information about a locus) and remove these from the dataset before mapping, using the <code>screen_for_duplicate_markers</code> function:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">screened_data4 &lt;-<span class="st"> </span><span class="kw">screen_for_duplicate_markers</span>(<span class="dt">dosage_matrix =</span> screened_data3)</code></pre></div>
<p><img src="" /><!-- --></p>
<pre><code>## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%
## Marked and merged 3 duplicated markers</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">dim</span>(screened_data4$filtered_dosage_matrix) </code></pre></div>
<pre><code>## [1] 1412  200</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">filtered_data &lt;-<span class="st"> </span>screened_data4$filtered_dosage_matrix </code></pre></div>
<p>Before we start the mapping, it is a good idea to run the summary function again, this time on the <code>filtered_data</code> object, so we have a clear record of the numbers of markers going forward to the mapping stage.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">pq_screened_data &lt;-<span class="st"> </span><span class="kw">parental_quantities</span>(<span class="dt">dosage_matrix =</span> filtered_data)</code></pre></div>
<p><img src="" /><!-- --></p>
<pre><code>##              1x0   2x0   0x1   1x1   2x1   0x2   1x2   2x2   1x3
## ----------  ----  ----  ----  ----  ----  ----  ----  ----  ----
## frequency    175   111   194   301   183   101   198    88    61</code></pre>
<p><a name="sect8"></a></p>
</div>
</div>
<div id="simplex-x-nulliplex-markers-defining-chromosomes-and-homologues" class="section level2">
<h2>8. Simplex x nulliplex markers – defining chromosomes and homologues</h2>
<p>The first step in actually mapping the markers is to define linkage groups for the markers. We use the simplex x nulliplex markers to define the linkage groups (chromosomes) and following that, the homologues (4 per linkage group are expected in a tetraploid). We use the LOD scores to cluster markers, which increases dramatically for tight coupling-phase estimates.</p>
<p>To estimate the recombination frequencies, LOD and phasing between all marker pairs of P1 enter:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SN_SN_P1 &lt;-<span class="st"> </span><span class="kw">linkage</span>(<span class="dt">dosage_matrix =</span> filtered_data, 
                    <span class="dt">markertype1 =</span> <span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">0</span>),
                    <span class="dt">target_parent =</span> <span class="st">&quot;P1&quot;</span>,
                    <span class="dt">other_parent =</span> <span class="st">&quot;P2&quot;</span>,
                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                    <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>
)</code></pre></div>
<blockquote>
<h6 id="a-note-on-arguments-and-multi-core-processing">A note on arguments and multi-core processing</h6>
<p>There is much more to the <code>linkage</code> function. Checkout <code>?linkage</code> to see the arguments that can be passed to it. One important feature is that it can run in parallel on multiple cores. This can be very time-efficient if you are working with a large number of markers. The number of cores to use can be specified using the <code>ncores</code> argument. Before doing so, first check the number of cores on your system using <code>parallel::detectCores()</code>. Use at least one core less than you have available to prevent overcommitment.</p>
</blockquote>
<p>It is a good idea to first visualise the relationship between the recombination frequency and LOD score, which gives us some insight into the differences between the information content of the different phases, as well as highlighting possible problems in the data. This can be done by using the <code>plot_linkage_df</code> function:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot_linkage_df</span>(SN_SN_P1, <span class="dt">r_max =</span> <span class="fl">0.5</span>)</code></pre></div>
<p><img src="" /><!-- --></p>
<p>In some situations, it is possible to identify outlying datapoints on this plot, which do not fit the expected relationship. We have implemented a function to check this for the case of pairs of simplex x nulliplex markers, as these are the markers that we use in the subsequent steps of clustering and identifying chromosome and homologue linkage groups. If you are unsure, you can run the <code>SNSN_LOD_deviation</code> function anyway to see whether your data fits the expected pattern:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P1deviations &lt;-<span class="st"> </span><span class="kw">SNSN_LOD_deviations</span>(<span class="dt">linkage_df =</span> SN_SN_P1,
                                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                    <span class="dt">N =</span> <span class="kw">ncol</span>(filtered_data) -<span class="st"> </span><span class="dv">2</span>, <span class="co">#The F1 population size</span>
                                    <span class="dt">alpha =</span> <span class="kw">c</span>(<span class="fl">0.05</span>,<span class="fl">0.2</span>),
                                    <span class="dt">plot_expected =</span> <span class="ot">TRUE</span>,
                                    <span class="dt">phase=</span><span class="st">&quot;coupling&quot;</span>)</code></pre></div>
<p><img src="" /><!-- --></p>
<p>Here we see there are a few “outliers” in the coupling pairwise data (identified by the stars) from the expected relationship between r and LOD (the bounds of which are shown by the dotted lines - these can be widened as necessary using the <code>alpha</code> argument for the upper and lower tolerances around the “true” line, usually in a trial and error approach). In our example dataset, these highlighted points are hardly outliers and we can proceed. In other datasets, you may want to remove these starred pairwise estimates before continuing with marker clustering, as it is possible they may cause unrelated homologue clusters to clump together and prevent a clear clustering structure from being defined. To do so, remove them from the linkage data frame as follows:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SN_SN_P1<span class="fl">.1</span> &lt;-<span class="st"> </span>SN_SN_P1[SN_SN_P1$phase ==<span class="st"> &quot;coupling&quot;</span>,][-<span class="kw">which</span>(P1deviations &gt;<span class="st"> </span><span class="fl">0.2</span>),]</code></pre></div>
<p>We are now in a position to cluster the marker data into chromosomes using the <code>cluster_SN_markers</code> function. This function clusters markers over a range of LOD thresholds. In the example below this is over a range from LOD 3 to 10 in steps of 1 (<code>LOD_sequence = c(3:10)</code>). If the argument <code>plot_network</code> is set to <code>TRUE</code> this function will plot the network of the lowest LOD threshold (LOD 3 in this case), and overlays the clusters of the higher LOD scores. This setting is not recommended with large number of marker pairs. Your computer will most likely run out of memory.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P1_homologues &lt;-<span class="st"> </span><span class="kw">cluster_SN_markers</span>(<span class="dt">linkage_df =</span> SN_SN_P1, 
                                    <span class="dt">LOD_sequence =</span> <span class="kw">seq</span>(<span class="dv">3</span>, <span class="dv">10</span>, <span class="fl">0.5</span>), 
                                    <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                    <span class="dt">parentname =</span> <span class="st">&quot;P1&quot;</span>,
                                    <span class="dt">plot_network =</span> <span class="ot">FALSE</span>,
                                    <span class="dt">plot_clust_size =</span> <span class="ot">FALSE</span>)</code></pre></div>
<pre><code>## Total number of edges: 894</code></pre>
<p><img src="" /><!-- --><img src="" /><!-- --></p>
<p>The function <code>cluster_SN_markers</code> returns a list of cluster data frames, one data frame for each LOD threshold. The first plot shows two things - on the normal y-axis, the number of clusters over a range of LOD thresholds (x-axis) - and on the right-hand y-axis, the number of markers that do not form part of any cluster, which is shown by the blue dashed line. In the <code>cluster_SN_markers</code> function, we specify the ploidy (4 for a tetraploid) and the expected number of chromosomes (LG_number = 5) and therefore there are 4 x 5 = 20 homologue clusters expected. If we examine the second plot, we see how these clusters stay together or fragment as the linkage stringency (LOD score) increases. Here, we see clearly that at LOD 3.5, there are 5 clusters which all eventually divide into 4 sub-clusters each. This is the ideal scenario - we have identified both the chromosomal linkage groups and the homologues using only the simplex x nulliplex marker data. Note that it is the coupling-phase estimates which define the homologue clusters and the repulsion-phase estimates which provide the associations between these coupling-clusters (homologues). In practice, things might not always be so orderly. In noisy datasets with higher numbers of genotyping errors, or at higher ploidy levels, the repulsion-phase linkages between homologues will be lost among false positive linkages between un-linked markers. This means that all markers will likely clump together at lower LOD scores.</p>
<p>You can manually check how many clusters there are at a specific LOD score, and how many markers are in a cluster:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P1_hom_LOD3<span class="fl">.5</span> &lt;-<span class="st"> </span>P1_homologues[[<span class="st">&quot;3.5&quot;</span>]]
t &lt;-<span class="st"> </span><span class="kw">table</span>(P1_hom_LOD3<span class="fl">.5</span>$cluster)
<span class="kw">print</span>(<span class="kw">paste</span>(<span class="st">&quot;Number of clusters:&quot;</span>,<span class="kw">length</span>(t)))</code></pre></div>
<pre><code>## [1] &quot;Number of clusters: 5&quot;</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">t[<span class="kw">order</span>(<span class="kw">as.numeric</span>(<span class="kw">names</span>(t)))]</code></pre></div>
<pre><code>## 
##  1  2  3  4  5 
## 37 36 29 31 41</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P1_hom_LOD5 &lt;-<span class="st"> </span>P1_homologues[[<span class="st">&quot;5&quot;</span>]]
t &lt;-<span class="st"> </span><span class="kw">table</span>(P1_hom_LOD5$cluster)
<span class="kw">print</span>(<span class="kw">paste</span>(<span class="st">&quot;Number of clusters:&quot;</span>,<span class="kw">length</span>(t)))</code></pre></div>
<pre><code>## [1] &quot;Number of clusters: 20&quot;</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">t[<span class="kw">order</span>(<span class="kw">as.numeric</span>(<span class="kw">names</span>(t)))]</code></pre></div>
<pre><code>## 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
## 13  9 10  4 12  7  9  8  3 13  6  7  4 10  6 11 17  6  6 12</code></pre>
<p>Based on the plots (or the output as shown) you should be able to decide which LOD threshold best separates the data. In this example, at a LOD of 3.5 we identify chromosomes and anywhere between LOD 4.5 and LOD 7 we have identified our homologues (4 x 5 = 20 clusters). In this parent (P1) we have therefore solved the marker clustering problem without any difficulty. It remains to gather this into a structure that can be used at later stages. To do this, we need two LOD values which define chromosomal groupings (here, <code>LOD_chm</code> = 3.5) and homologue groupings (here, <code>LOD_hom</code> = 5). We can then use the <code>define_LG_structure</code> function with these arguments as follows (note the argument LG_number refers to the expected number of chromosomes for this species):</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">LGHomDf_P1 &lt;-<span class="st"> </span><span class="kw">define_LG_structure</span>(<span class="dt">cluster_list =</span> P1_homologues, 
                                  <span class="dt">LOD_chm =</span> <span class="fl">3.5</span>, 
                                  <span class="dt">LOD_hom =</span> <span class="dv">5</span>, 
                                  <span class="dt">LG_number =</span> <span class="dv">5</span>)</code></pre></div>
<pre><code>## 
## ####SxN Marker(s) lost in clustering step at LOD 5:
## 
## |_           |
## |:-----------|
## |Ls_rs099908 |</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">head</span>(LGHomDf_P1)</code></pre></div>
<pre><code>##     SxN_Marker LG homologue
## 1  Ac_ns035629  1         1
## 3  Ac_ns045766  1         1
## 7  Ac_ns096474  1         1
## 22 Ac_ts072686  1         1
## 24 Ac_ws013795  1         1
## 42 Ap_rs075640  1         1</code></pre>
<p>However, you may find that it is not possible to identify chromosomal groups as simply as this. We therefore recommend that homologue groups be first identified, and then another marker type be used to provide associations between these clusters, and hence help form chromosomal linkage groups. This can be done as follows:</p>
<ul>
<li><p>Generate homologue clusters using the <code>cluster_SN_markers</code> function at a high LOD.</p></li>
<li><p>Use the Duplex x Nulliplex (or Simplex x Simplex) marker information to put the clusters back together into chromosomal groups.</p></li>
</ul>
<p>It makes sense to use only SxN marker pairs in coupling phase for the homologue clustering (since these associations define the homologues). Note that the most likely phase per marker pair (coupling or repulsion) was already estimated by our previous call to the <code>linkage</code> function.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SN_SN_P1_coupl &lt;-<span class="st"> </span>SN_SN_P1[SN_SN_P1$phase ==<span class="st"> &quot;coupling&quot;</span>,] <span class="co"># select only markerpairs in coupling</span>

P1_homologues_1 &lt;-<span class="st"> </span><span class="kw">cluster_SN_markers</span>(<span class="dt">linkage_df =</span> SN_SN_P1_coupl, 
                                    <span class="dt">LOD_sequence =</span> <span class="kw">c</span>(<span class="dv">3</span>:<span class="dv">12</span>), 
                                    <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                    <span class="dt">parentname =</span> <span class="st">&quot;P1&quot;</span>,
                                    <span class="dt">plot_network =</span> <span class="ot">FALSE</span>,
                                    <span class="dt">plot_clust_size =</span> <span class="ot">FALSE</span>)</code></pre></div>
<pre><code>## Total number of edges: 682</code></pre>
<p><img src="" /><!-- --><img src="" /><!-- --></p>
<p>In the dot plot above, we see that all chromosomal associations have now been removed with the repulsion information.</p>
<blockquote>
<h6 id="a-note-on-nested-lists">A note on nested lists</h6>
<p>In the <code>polymapR</code> package, output is often given as a nested list. Some basic information on lists is available on e.g. <a href="http://www.r-tutor.com/r-introduction/list">r-tutor</a>. We use nested lists because we are often dealing with hierarchical data: marker data &gt; homologue &gt; linkage group &gt; organism. Some basic R knowledge is needed to retrieve, visualize or filter data from such a nested list. However, if you would like to visualize it without R (e.g. in your favourite spreadsheet program), you can choose to export a nested list using <code>write_nested_list</code>. It will keep its hierarchical structure by writing it to directories that have the same structure as the the nested list. Here’s an example: <code>write_nested_list(list = P1_homologues, directory = &quot;P1_homologues&quot;)</code></p>
</blockquote>
<p>Next, we use linkage to the DxN markers to provide bridging linkages between homologue clusters, using the <code>bridgeHomologues</code> function. The first time we run this function, we use a LOD threshold of 4 for evidence of linkage between a SxN and a DxN marker. However, there may also be some false positive linkages at this level. It might therefore be necessary to run the function again at a higher LOD (say LOD 8) to avoid such false positives. In doing so, we might lose true linkage information. Therefore, it’s a good idea to run the function at LOD 4 first, then increase to a higher LOD if necessary, and build up a picture of which clusters form part of the same chromosome. First calculate linkage between SxN and DxN markers:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SN_DN_P1 &lt;-<span class="st"> </span><span class="kw">linkage</span>(<span class="dt">dosage_matrix =</span> filtered_data, 
                    <span class="dt">markertype1 =</span> <span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">0</span>),
                    <span class="dt">markertype2 =</span> <span class="kw">c</span>(<span class="dv">2</span>,<span class="dv">0</span>),
                    <span class="dt">target_parent =</span> <span class="st">&quot;P1&quot;</span>,
                    <span class="dt">other_parent =</span> <span class="st">&quot;P2&quot;</span>,
                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                    <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>)</code></pre></div>
<p>Then, use <code>bridgeHomologues</code> to use these SxN - DxN linkages to find associations (“bridges”) between homologue clusters. Note that we must provide a particular division of the Simplex x Nulliplex data - in this example, we use the division at LOD 6 which gave us 20 clusters (our 20 putative homologues):</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">LGHomDf_P1_1 &lt;-<span class="st"> </span><span class="kw">bridgeHomologues</span>(<span class="dt">cluster_stack =</span> P1_homologues_1[[<span class="st">&quot;6&quot;</span>]], 
                               <span class="dt">linkage_df =</span> SN_DN_P1, 
                               <span class="dt">LOD_threshold =</span> <span class="dv">4</span>, 
                               <span class="dt">automatic_clustering =</span> <span class="ot">TRUE</span>, 
                               <span class="dt">LG_number =</span> <span class="dv">5</span>,
                               <span class="dt">parentname =</span> <span class="st">&quot;P1&quot;</span>)</code></pre></div>
<p><img src="" /><!-- --></p>
<p>This results in an ideal situation of 5 linkage groups each with exactly 4 homologues:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">table</span>(LGHomDf_P1_1$LG, LGHomDf_P1_1$homologue)</code></pre></div>
<pre><code>##    
##      1  2  3  4
##   1 17 12  6  6
##   2 12  7  9  8
##   3  3 13  5  7
##   4  4 10  6 11
##   5 13  9  4 10</code></pre>
<p>However, clustering of P2 is a bit less ideal. The following redoes everything for P2:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SN_SN_P2 &lt;-<span class="st"> </span><span class="kw">linkage</span>(<span class="dt">dosage_matrix =</span> filtered_data, 
                    <span class="dt">markertype1 =</span> <span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">0</span>),
                    <span class="dt">target_parent =</span> <span class="st">&quot;P2&quot;</span>,
                    <span class="dt">other_parent =</span> <span class="st">&quot;P1&quot;</span>,
                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                    <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>
)</code></pre></div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SN_SN_P2_coupl &lt;-<span class="st"> </span>SN_SN_P2[SN_SN_P2$phase ==<span class="st"> &quot;coupling&quot;</span>,] <span class="co"># get only markerpairs in coupling</span>

P2_homologues &lt;-<span class="st"> </span><span class="kw">cluster_SN_markers</span>(<span class="dt">linkage_df =</span> SN_SN_P2_coupl, 
                                    <span class="dt">LOD_sequence =</span> <span class="kw">c</span>(<span class="dv">3</span>:<span class="dv">12</span>), 
                                    <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                    <span class="dt">parentname =</span> <span class="st">&quot;P2&quot;</span>,
                                    <span class="dt">plot_network =</span> <span class="ot">FALSE</span>,
                                    <span class="dt">plot_clust_size =</span> <span class="ot">FALSE</span>)</code></pre></div>
<pre><code>## Total number of edges: 743</code></pre>
<p><img src="" /><!-- --><img src="" /><!-- --></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SN_DN_P2 &lt;-<span class="st"> </span><span class="kw">linkage</span>(<span class="dt">dosage_matrix =</span> filtered_data, 
                    <span class="dt">markertype1 =</span> <span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">0</span>),
                    <span class="dt">markertype2 =</span> <span class="kw">c</span>(<span class="dv">2</span>,<span class="dv">0</span>),
                    <span class="dt">target_parent =</span> <span class="st">&quot;P2&quot;</span>,
                    <span class="dt">other_parent =</span> <span class="st">&quot;P1&quot;</span>,
                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                    <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>)</code></pre></div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">LGHomDf_P2 &lt;-<span class="st"> </span><span class="kw">bridgeHomologues</span>(<span class="dt">cluster_stack =</span> P2_homologues[[<span class="st">&quot;6&quot;</span>]], 
                               <span class="dt">linkage_df =</span> SN_DN_P2, 
                               <span class="dt">LOD_threshold =</span> <span class="dv">4</span>, 
                               <span class="dt">automatic_clustering =</span> <span class="ot">TRUE</span>, 
                               <span class="dt">LG_number =</span> <span class="dv">5</span>,
                               <span class="dt">parentname =</span> <span class="st">&quot;P2&quot;</span>)</code></pre></div>
<p><img src="" /><!-- --></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">table</span>(LGHomDf_P2$LG,LGHomDf_P2$homologue)</code></pre></div>
<pre><code>##    
##      1  2  3  4  5
##   1 11  9  3  5  0
##   2 13 11 11  8  2
##   3 10 10  6  9  0
##   4 13 13  7 11  0
##   5  9  9 13 11  0</code></pre>
<p>Linkage group 2 possesses five “homologues” instead of four - it is therefore likely we have a number of fragmented homologues. However, we don’t know which clusters are actual homologues, and which are fragments. We can try to determine this using two methods.</p>
<div id="using-the-function-overviewsnlinks" class="section level3">
<h3>8.1 Using the function <code>overviewSNlinks</code></h3>
<p>Run <code>overviewSNlinks</code> as follows:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">overviewSNlinks</span>(<span class="dt">linkage_df =</span> SN_SN_P2,
                <span class="dt">LG_hom_stack =</span> LGHomDf_P2,
                <span class="dt">LG_number =</span> <span class="dv">2</span>,
                <span class="dt">LOD_threshold =</span> <span class="dv">2</span>)</code></pre></div>
<p><img src="" /><!-- --></p>
<p>Apparently, homologue 4 and 5 have markers in coupling, whereas all other homologue combinations have only markers in repulsion. We therefore merge homologue 4 and 5 (this is specified in the argument <code>mergeList</code> - use <code>?merge_homologues</code> for more information):</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">LGHomDf_P2_1 &lt;-<span class="st"> </span><span class="kw">merge_homologues</span>(<span class="dt">LG_hom_stack =</span> LGHomDf_P2,
                                 <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                 <span class="dt">linkage_group =</span> <span class="dv">2</span>,
                                 <span class="dt">mergeList =</span> <span class="kw">list</span>(<span class="kw">c</span>(<span class="dv">4</span>,<span class="dv">5</span>)))</code></pre></div>
<pre><code>## ####Number of markers per homologue:
## 
##        1    2    3    4    5
## ---  ---  ---  ---  ---  ---
## 1     11   13   10   13    9
## 2      9   11   10   13    9
## 3      3   11    6    7   13
## 4      5   10    9   11   11</code></pre>
</div>
<div id="using-the-function-cluster_per_lg" class="section level3">
<h3>8.2 Using the function <code>cluster_per_LG</code></h3>
<p><code>cluster_per_LG</code> does not rely on markers in repulsion, but uses variation in LOD thresholds between markers in coupling. Run it as follows:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">cluster_per_LG</span>(<span class="dt">LG =</span> <span class="dv">2</span>,
               <span class="dt">linkage_df =</span> SN_SN_P2[SN_SN_P2$phase ==<span class="st"> &quot;coupling&quot;</span>,], 
               <span class="dt">LG_hom_stack =</span> LGHomDf_P2, 
               <span class="dt">LOD_sequence =</span> <span class="kw">c</span>(<span class="dv">3</span>:<span class="dv">10</span>), <span class="co"># The first element is used for network layout</span>
               <span class="dt">modify_LG_hom_stack =</span> <span class="ot">FALSE</span>, 
               <span class="dt">network.layout =</span> <span class="st">&quot;stacked&quot;</span>,
               <span class="dt">nclust_out =</span> <span class="dv">4</span>,
               <span class="dt">label.offset=</span><span class="fl">1.2</span>)</code></pre></div>
<pre><code>## Total number of edges: 185</code></pre>
<p><img src="" /><!-- --></p>
<p>The network layout is based on the first LOD score you provide in <code>LOD_sequence</code> (which is LOD = 3 in this example). The darker the background, the more consistent the clustering is over the range of LOD scores. Change the <code>LOD_sequence</code> if the network is not as required (4 clusters in the case of a tetraploid) and re-run. If the network looks as you would like to, use the <code>cluster_per_LG</code> function but with a <code>LOD_sequence</code> of length 1 (the first LOD score, so again LOD = 3 in this example), and specify <code>modify_LG_hom_stack = TRUE</code> to make a modification to the cluster numbering:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">LGHomDf_P2_1 &lt;-<span class="st"> </span><span class="kw">cluster_per_LG</span>(<span class="dt">LG =</span> <span class="dv">2</span>, 
                               <span class="dt">linkage_df =</span> SN_SN_P2[SN_SN_P2$phase ==<span class="st"> &quot;coupling&quot;</span>,], 
                               <span class="dt">LG_hom_stack =</span> LGHomDf_P2, 
                               <span class="dt">LOD_sequence =</span> <span class="dv">3</span>, 
                               <span class="dt">modify_LG_hom_stack =</span> <span class="ot">TRUE</span>, 
                               <span class="dt">network.layout =</span> <span class="st">&quot;n&quot;</span>,
                               <span class="dt">nclust_out =</span> <span class="dv">4</span>)</code></pre></div>
<pre><code>## Total number of edges: 185</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">table</span>(LGHomDf_P2_1$homologue, LGHomDf_P2_1$LG)</code></pre></div>
<pre><code>##    
##      1  2  3  4  5
##   1 11 10 10 13  9
##   2  9 13 10 13  9
##   3  3 11  6  7 13
##   4  5 11  9 11 11</code></pre>
</div>
<div id="cross-ploidy-populations-e.g.tetraploid-x-diploid-to-give-triploid-f1" class="section level3">
<h3>8.3 Cross-ploidy populations (e.g. tetraploid x diploid to give triploid F1)</h3>
<p><code>polymapR</code> is also able to map triploid populations, and may be extended to other cross-ploidy populations if there is demand for that in future. However, because of differences in parental ploidy levels, we have implemented a mapping approach that assumes that the parent of higher ploidy level is parent 1. Usually in <code>polymapR</code> a cross is understood to be between “maternal parent” x “paternal parent”. However this is just a convenient convention, and parent numbering can be reversed if so desired. In the case of e.g. triploid populations, there is no flexibility anymore, parent 1 <em>must</em> be tetraploid and parent 2 <em>must</em> be diploid.</p>
<p>The next point to be aware of is that the segregation of markers in the diploid parent is assumed to follow disomic inheritance (there is no other option in a diploid), whereas the segregation of the tetraploid parent is assumed to follow tetrasomic inheritance. Mixtures with preferential pairing in the tetraploid parent are currently not implemented for triploid populations.</p>
<p>Because of this disomic behaviour in the diploid parent, we are unable to use the normal clustering function (<code>cluster_SN_markers</code>) for identifying the homologues. Proceed as normal with the steps described above for parent 1. However in order to proceed correctly for parent 2 (diploid), you will need to provide the markertype as <code>c(0,1)</code> and not <code>c(1,0)</code>. This is easy to understand since a 1x0 marker in a triploid population is <strong>Aaaa x aa</strong> while a 0x1 marker is <strong>aaaa x Aa</strong>. These are clearly different marker types and we cannot reverse them. To access the correct linkage functions for the diploid parent, therefore, you will need to think in terms of 1x0 markers segregating from parent 1, and 0x1 markers segregating from parent 2.</p>
<p>If you are confused about what marker types are possible in a triploid population, run the following, which gives all possible marker combinations (and the segregation ratios in the F1):</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">get</span>(<span class="st">&quot;seg_p3_random&quot;</span>,<span class="dt">envir=</span><span class="kw">getNamespace</span>(<span class="st">&quot;polymapR&quot;</span>))</code></pre></div>
<pre><code>##    dosage1 dosage2 f1_0 f1_1 f1_2 f1_3
## 1        0       0    1    0    0    0
## 2        0       1    1    1    0    0
## 3        0       2    0    1    0    0
## 4        1       0    1    1    0    0
## 5        1       1    1    2    1    0
## 6        1       2    0    1    1    0
## 7        2       0    1    4    1    0
## 8        2       1    1    5    5    1
## 9        2       2    0    1    4    1
## 10       3       0    0    1    1    0
## 11       3       1    0    1    2    1
## 12       3       2    0    0    1    1
## 13       4       0    0    0    1    0
## 14       4       1    0    0    1    1
## 15       4       2    0    0    0    1</code></pre>
<p>To run the preliminary linkage analysis of 0x1 markers in parent 2, we also need to specify the ploidy of parent 2, in this case <code>ploidy2 = 2</code>. Note that the target parent is P2 as this is the parent we are now interested in, but <code>markertype1 = c(0,1)</code>. For this we need a special dataset with triploid data, a sample of which is available from the <code>polymapR</code> package also:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">data</span>(<span class="st">&quot;TRI_dosages&quot;</span>)

SN_SN_P2.tri &lt;-<span class="st"> </span><span class="kw">linkage</span>(<span class="dt">dosage_matrix =</span> TRI_dosages, 
                    <span class="dt">markertype1 =</span> <span class="kw">c</span>(<span class="dv">0</span>,<span class="dv">1</span>),
                    <span class="dt">target_parent =</span> <span class="st">&quot;P2&quot;</span>, <span class="co">#Note the target parent is P2</span>
                    <span class="dt">other_parent =</span> <span class="st">&quot;P1&quot;</span>,
                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                    <span class="dt">ploidy2 =</span> <span class="dv">2</span>,
                    <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>
)</code></pre></div>
<p>Have a look at the plot of r versus LOD score - here you should notice that both the coupling and repulsion phases are equally informative now and cannot therefore be separated by LOD score as before:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot_linkage_df</span>(SN_SN_P2.tri)</code></pre></div>
<p><img src="" /><!-- --></p>
<p>If we try to cluster using the old function, we get the following result - no separation of the homologues, but the chromosomal linkage groups are identified:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P2_homologues.tri &lt;-<span class="st"> </span><span class="kw">cluster_SN_markers</span>(<span class="dt">linkage_df =</span> SN_SN_P2.tri, 
                                    <span class="dt">LOD_sequence =</span> <span class="kw">seq</span>(<span class="dv">3</span>, <span class="dv">10</span>, <span class="dv">1</span>), 
                                    <span class="dt">LG_number =</span> <span class="dv">3</span>,
                                    <span class="dt">ploidy =</span> <span class="dv">2</span>, <span class="co">#because P2 is diploid..</span>
                                    <span class="dt">parentname =</span> <span class="st">&quot;P2&quot;</span>,
                                    <span class="dt">plot_network =</span> <span class="ot">FALSE</span>,
                                    <span class="dt">plot_clust_size =</span> <span class="ot">FALSE</span>) </code></pre></div>
<pre><code>## Total number of edges: 345</code></pre>
<p><img src="" /><!-- --></p>
<pre><code>## Please note: the number of clusters is smaller than number of expected homologues.
##             Not all expected homologues are represented.</code></pre>
<p><img src="" /><!-- --></p>
<p>To separate the 0x1 markers into homologues, we use the phase information from the linkage analysis directly, using the <code>phase_SN_diploid</code> function. To learn more about how this function works, input <code>?phase_SN_diploid</code> first, to get an idea of the arguments used.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">LGHomDf_P2.tri &lt;-<span class="st"> </span><span class="kw">phase_SN_diploid</span>(<span class="dt">linkage_df =</span> SN_SN_P2.tri,
                                   <span class="dt">cluster_list =</span> P2_homologues.tri,
                                   <span class="dt">LOD_chm =</span> <span class="dv">4</span>, <span class="co">#LOD at which chromosomes are identified</span>
                                   <span class="dt">LG_number =</span> <span class="dv">3</span>) <span class="co">#number of linkage groups</span></code></pre></div>
<pre><code>## Total number of edges: 161
## Complete phase assignment possible using only coupling information at LOD 4</code></pre>
<p>The rest of the analysis should be pretty much the same as for tetraploids or hexaploids, although you should check each time you use a new function whether that function has a <code>ploidy2</code> argument, and if it does, <strong>use it!</strong></p>
<p><a name="sect9"></a></p>
</div>
</div>
<div id="assigning-sxs-and-dxn-markers-and-consensus-linkage-group-lg-names" class="section level2">
<h2>9. Assigning SxS and DxN markers and consensus linkage group (LG) names</h2>
<p>Assigning simplex x simplex markers to a homologue and linkage group is done by calculating linkage between SxS and SxN markers and after that assigning them based on this linkage using <code>assign_linkage_group</code>:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SN_SS_P1 &lt;-<span class="st"> </span><span class="kw">linkage</span>(<span class="dt">dosage_matrix =</span> filtered_data, 
                    <span class="dt">markertype1 =</span> <span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">0</span>),
                    <span class="dt">markertype2 =</span> <span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">1</span>),
                    <span class="dt">target_parent =</span> <span class="st">&quot;P1&quot;</span>,
                    <span class="dt">other_parent =</span> <span class="st">&quot;P2&quot;</span>,
                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                    <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>)</code></pre></div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P1_SxS_Assigned &lt;-<span class="st"> </span><span class="kw">assign_linkage_group</span>(<span class="dt">linkage_df =</span> SN_SS_P1,
                                        <span class="dt">LG_hom_stack =</span> LGHomDf_P1,
                                        <span class="dt">SN_colname =</span> <span class="st">&quot;marker_a&quot;</span>,
                                        <span class="dt">unassigned_marker_name =</span> <span class="st">&quot;marker_b&quot;</span>,
                                        <span class="dt">phase_considered =</span> <span class="st">&quot;coupling&quot;</span>,
                                        <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                        <span class="dt">LOD_threshold =</span> <span class="dv">3</span>,
                                        <span class="dt">ploidy =</span> <span class="dv">4</span>)</code></pre></div>
<pre><code>## 
## ####Marker(s) showing ambiguous linkage to more than one LG:
## 
## |_           |
## |:-----------|
## |Ap_ts069042 |
## 
##  In total, 296 out of 301 markers were assigned.
## 
## ####Marker(s) not assigned:
## 
## _             _             _             _           
## ------------  ------------  ------------  ------------
## Ac_ws033297   Ap_ts089267   St_ns048176   St_ts030296 
## Zm_ws010615</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">head</span>(P1_SxS_Assigned)</code></pre></div>
<pre><code>##             Assigned_LG LG1 LG2 LG3 LG4 LG5 Hom1 Hom2 Hom3 Hom4
## Ac_ns002135           4   0   0   0   7   0    0    7    0    0
## Ac_ns002510           2   0   7   0   0   0    0    0    0    7
## Ac_ns024650           2   0  11   0   0   0   11    0    0    0
## Ac_ns028519           5   0   0   0   0   5    0    5    0    0
## Ac_ns028533           5   0   0   0   0   3    0    0    3    0
## Ac_ns029178           3   0   0  12   0   0    0   12    0    0
##             Assigned_hom1 Assigned_hom2 Assigned_hom3 Assigned_hom4
## Ac_ns002135             2            NA            NA            NA
## Ac_ns002510             4            NA            NA            NA
## Ac_ns024650             1            NA            NA            NA
## Ac_ns028519             2            NA            NA            NA
## Ac_ns028533             3            NA            NA            NA
## Ac_ns029178             2            NA            NA            NA</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">SN_SS_P2 &lt;-<span class="st"> </span><span class="kw">linkage</span>(<span class="dt">dosage_matrix =</span> filtered_data, 
                    <span class="dt">markertype1 =</span> <span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">0</span>),
                    <span class="dt">markertype2 =</span> <span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">1</span>),
                    <span class="dt">target_parent =</span> <span class="st">&quot;P2&quot;</span>,
                    <span class="dt">other_parent =</span> <span class="st">&quot;P1&quot;</span>,
                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                    <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>)</code></pre></div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P2_SxS_Assigned &lt;-<span class="st"> </span><span class="kw">assign_linkage_group</span>(<span class="dt">linkage_df =</span> SN_SS_P2,
                                        <span class="dt">LG_hom_stack =</span> LGHomDf_P2_1,
                                        <span class="dt">SN_colname =</span> <span class="st">&quot;marker_a&quot;</span>,
                                        <span class="dt">unassigned_marker_name =</span> <span class="st">&quot;marker_b&quot;</span>,
                                        <span class="dt">phase_considered =</span> <span class="st">&quot;coupling&quot;</span>,
                                        <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                        <span class="dt">LOD_threshold =</span> <span class="dv">3</span>,
                                        <span class="dt">ploidy =</span> <span class="dv">4</span>)</code></pre></div>
<pre><code>## 
##  In total, 300 out of 300 markers were assigned.</code></pre>
<p>As simplex x simplex markers are present in both parents, we can define which linkage groups correspond with each other between parents. After this, one of the linkage groups of the parents should be renamed and the SxS markers assigned again according to the new linkage group names:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">LGHomDf_P2_2 &lt;-<span class="st"> </span><span class="kw">consensus_LG_names</span>(<span class="dt">modify_LG =</span> LGHomDf_P2_1, 
                                   <span class="dt">template_SxS =</span> P1_SxS_Assigned, 
                                   <span class="dt">modify_SxS =</span> P2_SxS_Assigned)</code></pre></div>
<pre><code>## 
## ####Original LG names
## 
##        1    2    3    4    5
## ---  ---  ---  ---  ---  ---
## 1      0    0    0    0   47
## 2     64    0    0    0    0
## 3      0   58    0    0    0
## 4      0    0    0   58    0
## 5      0    0   68    0    0
## 
## ####Modified LG names
## 
##        1    2    3    4    5
## ---  ---  ---  ---  ---  ---
## 1     47    0    0    0    0
## 2      0   64    0    0    0
## 3      0    0   58    0    0
## 4      0    0    0   58    0
## 5      0    0    0    0   68</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P2_SxS_Assigned &lt;-<span class="st"> </span><span class="kw">assign_linkage_group</span>(<span class="dt">linkage_df =</span> SN_SS_P2,
                                        <span class="dt">LG_hom_stack =</span> LGHomDf_P2_2,
                                        <span class="dt">SN_colname =</span> <span class="st">&quot;marker_a&quot;</span>,
                                        <span class="dt">unassigned_marker_name =</span> <span class="st">&quot;marker_b&quot;</span>,
                                        <span class="dt">phase_considered =</span> <span class="st">&quot;coupling&quot;</span>,
                                        <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                        <span class="dt">LOD_threshold =</span> <span class="dv">3</span>,
                                        <span class="dt">ploidy =</span> <span class="dv">4</span>)</code></pre></div>
<pre><code>## 
##  In total, 300 out of 300 markers were assigned.</code></pre>
<p>Since we now have a consistent linkage group numbering, we can also assign the DxN markers:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P1_DxN_Assigned &lt;-<span class="st"> </span><span class="kw">assign_linkage_group</span>(<span class="dt">linkage_df =</span> SN_DN_P1,
                                        <span class="dt">LG_hom_stack =</span> LGHomDf_P1,
                                        <span class="dt">SN_colname =</span> <span class="st">&quot;marker_a&quot;</span>,
                                        <span class="dt">unassigned_marker_name =</span> <span class="st">&quot;marker_b&quot;</span>,
                                        <span class="dt">phase_considered =</span> <span class="st">&quot;coupling&quot;</span>,
                                        <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                        <span class="dt">LOD_threshold =</span> <span class="dv">3</span>,
                                        <span class="dt">ploidy =</span> <span class="dv">4</span>)</code></pre></div>
<pre><code>## 
##  In total, 111 out of 111 markers were assigned.</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P2_DxN_Assigned &lt;-<span class="st"> </span><span class="kw">assign_linkage_group</span>(<span class="dt">linkage_df =</span> SN_DN_P2,
                                        <span class="dt">LG_hom_stack =</span> LGHomDf_P2_2,
                                        <span class="dt">SN_colname =</span> <span class="st">&quot;marker_a&quot;</span>,
                                        <span class="dt">unassigned_marker_name =</span> <span class="st">&quot;marker_b&quot;</span>,
                                        <span class="dt">phase_considered =</span> <span class="st">&quot;coupling&quot;</span>,
                                        <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                        <span class="dt">LOD_threshold =</span> <span class="dv">3</span>,
                                        <span class="dt">ploidy =</span> <span class="dv">4</span>)</code></pre></div>
<pre><code>## 
##  In total, 101 out of 101 markers were assigned.</code></pre>
<p><a name="sect10"></a></p>
</div>
<div id="assign-all-other-markertypes" class="section level2">
<h2>10. Assign all other markertypes</h2>
<p>Now that we have a backbone of SxN markers with consistent linkage group names, it is time to assign all other marker types to a linkage group and homologue using their linkage to SxN markers. Since we already did this for DxN and SxS markers, there is no need to re-do this work for these marker types. The function <code>homologue_lg_assignment</code> finds all markertypes that have not been assigned yet, does the linkage analysis and assigns them to a linkage group and homologue. As linkage is calculated for a lot of marker combinations, this might take a while.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">marker_assignments_P1 &lt;-<span class="st"> </span><span class="kw">homologue_lg_assignment</span>(<span class="dt">dosage_matrix =</span> filtered_data,
                                                 <span class="dt">assigned_list =</span> <span class="kw">list</span>(P1_SxS_Assigned, 
                                                                      P1_DxN_Assigned),
                                                 <span class="dt">assigned_markertypes =</span> <span class="kw">list</span>(<span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">1</span>), <span class="kw">c</span>(<span class="dv">2</span>,<span class="dv">0</span>)),
                                                 <span class="dt">LG_hom_stack =</span> LGHomDf_P1,
                                                 <span class="dt">target_parent =</span> <span class="st">&quot;P1&quot;</span>,
                                                 <span class="dt">other_parent =</span> <span class="st">&quot;P2&quot;</span>,
                                                 <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                                 <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>,
                                                 <span class="dt">convert_palindrome_markers =</span> <span class="ot">FALSE</span>,
                                                 <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                                 <span class="dt">LOD_threshold =</span> <span class="dv">3</span>,
                                                 <span class="dt">write_intermediate_files =</span> <span class="ot">FALSE</span>
)</code></pre></div>
<p>Also do this for P2:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">marker_assignments_P2 &lt;-<span class="st"> </span><span class="kw">homologue_lg_assignment</span>(<span class="dt">dosage_matrix =</span> filtered_data,
                                                 <span class="dt">assigned_list =</span> <span class="kw">list</span>(P2_SxS_Assigned, 
                                                                      P2_DxN_Assigned),
                                                 <span class="dt">assigned_markertypes =</span> <span class="kw">list</span>(<span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">1</span>), <span class="kw">c</span>(<span class="dv">2</span>,<span class="dv">0</span>)),
                                                 <span class="dt">LG_hom_stack =</span> LGHomDf_P2_2,
                                                 <span class="dt">target_parent =</span> <span class="st">&quot;P2&quot;</span>,
                                                 <span class="dt">other_parent =</span> <span class="st">&quot;P1&quot;</span>,
                                                 <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                                 <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>,
                                                 <span class="dt">convert_palindrome_markers =</span> <span class="ot">TRUE</span>,
                                                 <span class="dt">LG_number =</span> <span class="dv">5</span>,
                                                 <span class="dt">LOD_threshold =</span> <span class="dv">3</span>,
                                                 <span class="dt">write_intermediate_files =</span> <span class="ot">FALSE</span>
)</code></pre></div>
<p>Next, to make sure the marker linkage group assignment is consistent across parents we run the <code>check_marker_assignment</code> function, removing any bi-parental markers if they show linkage to different chromosomes (which suggests a problem with these markers):</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">marker_assignments &lt;-<span class="st"> </span><span class="kw">check_marker_assignment</span>(marker_assignments_P1,marker_assignments_P2)</code></pre></div>
<p><a name="sect11"></a></p>
</div>
<div id="finish-the-linkage-analysis" class="section level2">
<h2>11. Finish the linkage analysis</h2>
<p>Since all markers have been assigned to a linkage group, we now do the linkage calculations per linkage group of every marker type combination (so not only with SxN markers). This is done with the <code>finish_linkage_analysis</code> function:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">all_linkages_list_P1 &lt;-<span class="st"> </span><span class="kw">finish_linkage_analysis</span>(<span class="dt">marker_assignment =</span> marker_assignments$P1,
                                                <span class="dt">dosage_matrix =</span> filtered_data,
                                                <span class="dt">target_parent =</span> <span class="st">&quot;P1&quot;</span>,
                                                <span class="dt">other_parent =</span> <span class="st">&quot;P2&quot;</span>,
                                                <span class="dt">convert_palindrome_markers =</span> <span class="ot">FALSE</span>,
                                                <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                                <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>,
                                                <span class="dt">LG_number =</span> <span class="dv">5</span>) 

all_linkages_list_P2 &lt;-<span class="st"> </span><span class="kw">finish_linkage_analysis</span>(<span class="dt">marker_assignment =</span> marker_assignments$P2,
                                                <span class="dt">dosage_matrix =</span> filtered_data,
                                                <span class="dt">target_parent =</span> <span class="st">&quot;P2&quot;</span>,
                                                <span class="dt">other_parent =</span> <span class="st">&quot;P1&quot;</span>,
                                                <span class="dt">convert_palindrome_markers =</span> <span class="ot">TRUE</span>, <span class="co"># convert 3.1 markers</span>
                                                <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                                <span class="dt">pairing =</span> <span class="st">&quot;random&quot;</span>,
                                                <span class="dt">LG_number =</span> <span class="dv">5</span>)</code></pre></div>
<p>The output is returned in a list:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">str</span>(all_linkages_list_P1)</code></pre></div>
<p>Using the assignment data, we can also split this list in homologues using <code>split_linkage_info</code>:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">all_linkages_list_P1_split &lt;-<span class="st"> </span><span class="kw">split_linkage_info</span>(<span class="dt">all_linkages =</span> all_linkages_list_P1, 
                                                 <span class="dt">marker_assignment =</span> marker_assignments_P1, 
                                                 <span class="dt">ploidy =</span> <span class="dv">4</span>)

all_linkages_list_P2_split &lt;-<span class="st"> </span><span class="kw">split_linkage_info</span>(<span class="dt">all_linkages =</span> all_linkages_list_P2, 
                                                 <span class="dt">marker_assignment =</span> marker_assignments_P2, 
                                                 <span class="dt">ploidy =</span> <span class="dv">4</span>)</code></pre></div>
<p>This results in a nested list, which you can write to files:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">str</span>(all_linkages_list_P1_split)
<span class="kw">write_nested_list</span>(<span class="dt">nested_list =</span> all_linkages_list_P1, <span class="dt">directory =</span> <span class="st">&quot;linkages_P1&quot;</span>)</code></pre></div>
<p><a name="sect12"></a></p>
</div>
<div id="optional-marker-binning" class="section level2">
<h2>12. <em>Optional</em>: Marker binning</h2>
<p>If you are using JoinMap to do the ordering, the number of markers than can be used per homologue is limited. If the number of markers per homologue approaches or exceeds 100, we first need to perform some sort of binning strategy, putting co-segregating markers into “bins” and only using one of them as a representative marker for the others in the bin.</p>
<p>Marker binning can be done on one <code>linkage_df</code> separately, using the <code>marker_binning</code> function, or by using its wrapper <code>marker_binning_list</code> which performs binning across all linkage groups:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">all_linkages_list_P1.binned &lt;-<span class="st"> </span><span class="kw">marker_binning_list</span>(<span class="dt">dosage_matrix =</span> filtered_data, 
                                            <span class="dt">linkage_list =</span> all_linkages_list_P1_split,
                                            <span class="dt">target_parent =</span> <span class="st">&quot;P1&quot;</span>, 
                                            <span class="dt">other_parent =</span> <span class="st">&quot;P2&quot;</span>,
                                            <span class="dt">return_removed_marker_info =</span> <span class="ot">TRUE</span>)

all_linkages_list_P2.binned &lt;-<span class="st"> </span><span class="kw">marker_binning_list</span>(<span class="dt">dosage_matrix =</span> filtered_data, 
                                            <span class="dt">linkage_list =</span> all_linkages_list_P2_split,
                                            <span class="dt">target_parent =</span> <span class="st">&quot;P2&quot;</span>, 
                                            <span class="dt">other_parent =</span> <span class="st">&quot;P1&quot;</span>,
                                            <span class="dt">return_removed_marker_info =</span> <span class="ot">TRUE</span>)</code></pre></div>
<p><a name="sect13"></a></p>
</div>
<div id="optional-writing-to-a-.pwd-file" class="section level2">
<h2>13. <em>Optional</em>: Writing to a .pwd file</h2>
<p>If the marker ordering is to be done by the JoinMap GUI, the linkage data can be written to JoinMap compatible .pwd files using <code>write.pwd</code> or its wrapper <code>write_pwd_list</code>:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">write_pwd_list</span>(<span class="dt">linkages_list =</span> all_linkages_list_P1$linkage_list, 
               <span class="dt">target_parent =</span> <span class="st">&quot;P1&quot;</span>, 
               <span class="dt">binned =</span> T)

<span class="kw">write_pwd_list</span>(<span class="dt">linkages_list =</span> all_linkages_list_P2$linkage_list, 
               <span class="dt">target_parent =</span> <span class="st">&quot;P2&quot;</span>, 
               <span class="dt">binned =</span> T)</code></pre></div>
<p><a name="sect14"></a></p>
</div>
<div id="marker-ordering" class="section level2">
<h2>14. Marker ordering</h2>
<p>The basis of producing linkage maps has now been achieved, that is the pairwise recombination frequencies have been estimated and LOD scores for these estimates have been calculated for all markers. There are two options in the <code>polymapR</code> package for marker ordering: weighted linear regression (based on the JoinMap algorithm) or multi-dimensional scaling. The first is a reliable and widely used algorithm, but is extremely slow for higher marker numbers. The latter option has recently been developed by Katherine Preedy and Christine Hackett <span class="citation">[2]</span> and offers similar levels of accuracy as the weighted linear regression but at much higher speed. Our approach is to create genetic maps of the homologues as well as an integrated map across all hoologues and parents. As an example we use the multidimensional scaling approach, using the function <code>MDSMap_from_list</code>. This function applies the <code>estimate.map</code> function from the <code>MDSMap</code> package to a list of linkage dataframes. Below, we are using the default settings of <code>estimate.map</code>, however, they can be changed by supplying extra arguments.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">maplist_P1_LG1 &lt;-<span class="st"> </span><span class="kw">MDSMap_from_list</span>(all_linkages_list_P1_split$LG1, 
                               <span class="dt">write_to_file =</span> <span class="ot">FALSE</span>
                               )
maplist_P2_LG1 &lt;-<span class="st"> </span><span class="kw">MDSMap_from_list</span>(all_linkages_list_P2_split$LG1, 
                               <span class="dt">write_to_file =</span> <span class="ot">FALSE</span>)</code></pre></div>
<p>By calling the <code>MDSMap_from_list</code> function you automatically produce the input .txt files that <code>MDSMap</code> requires in the correct format. Note that unless you provide a different folder name (by specifying <code>mapdir</code>), the files will be overwritten each time you run the function. We generally recommend using the default mapping settings of the <code>estimate.map</code> function (so using the method of principal curves in 2 dimensions with LOD<sup>2</sup> as weights), but there are a number of options that can be specified which are described in the manual of that package (check out <a href="https://cran.R-project.org/package=MDSMap" class="uri">https://cran.R-project.org/package=MDSMap</a>). Output plots can be saved as .pdf files by specifying <code>write_to_file = TRUE</code>, in the same <code>mapdir</code> folder that contains map input files as well.</p>
</div>
<div id="creating-an-integrated-chromosomal-linkage-map" class="section level2">
<h2>14.1 Creating an integrated chromosomal linkage map</h2>
<p>In the previous example we created separate linkage maps per homologue. However, it is also possible to create an integrated chromosomal map using <code>MDSMap</code>, by simply supplying all the linkage information for the whole linkage group together. To create an integrated map of all linkage groups, we first combine the linkage information together and then run the <code>MDSMap_from_list</code> function which prepares the files and passes them to <code>MDSMap</code> for the mapping:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">linkages &lt;-<span class="st"> </span><span class="kw">list</span>()
for(lg in <span class="kw">names</span>(all_linkages_list_P1)){
  linkages[[lg]] &lt;-<span class="st"> </span><span class="kw">rbind</span>(all_linkages_list_P1[[lg]], all_linkages_list_P2[[lg]])
}

integrated.maplist &lt;-<span class="st"> </span><span class="kw">MDSMap_from_list</span>(linkages)</code></pre></div>
</div>
<div id="phasing-an-integrated-map" class="section level2">
<h2>14.2 Phasing an integrated map</h2>
<p>One final step that is useful is to generate phased linkage maps from the integrated linkage map and marker assignments. This provides information on the coverage of markers across the parental homologues. To phase the markers, we use the <code>create_phased_maplist</code> function as follows:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">phased.maplist &lt;-<span class="st"> </span><span class="kw">create_phased_maplist</span>(<span class="dt">maplist =</span> integrated.maplist,
                                        <span class="dt">dosage_matrix =</span> filtered_data,
                                        <span class="dt">N_linkages =</span> <span class="dv">5</span>,
                                        <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                        <span class="dt">marker_assignment.1 =</span> marker_assignments$P1,
                                        <span class="dt">marker_assignment.2 =</span> marker_assignments$P2)</code></pre></div>
<p>The <code>N_linkages</code> option specifies the minimum number of significant linkages of a marker to a chromosomal linkage group to be confident of its assignment. Note that this level of significance (e.g. LOD &gt; 3) was already defined in the <code>homologue_lg_assignment</code> function earlier. There are a number of other arguments with <code>create_phased_maplist</code>, including <code>original_coding</code> which produces a phased mapfile in the original, uncoverted format. This may have advantages for tracking marker alleles, although it has the disadvantage that the homologues appear to be more saturated with markers than they actually are!</p>
</div>
<div id="optional-independent-homologue-map-integration" class="section level2">
<h2>14.3 <em>Optional</em>: Independent homologue map integration</h2>
<p>In some situations it may be desirable to independently integrate homologue maps rather than directly generating integrated chromosomal maps as we did in the previous section. For example, there may be cases where the rate of recombination is different between parents and you would like to compare integrated chromosomal maps using different methods.</p>
<p>The <code>LPMerge</code> package <span class="citation">[3]</span> in <code>R</code> is able to integrate homologue maps in a polyploid but has a number of limitations, one of them being that it is unable to recognise when the maps to be integrated are not orientated correctly. A genetic position of 0 cM merely means we are at the outermost marker position in a linkage group – it does not correspond to a “0” position in any real sense. We therefore have to first test whether the homologue maps are consistently ordered before proceeding with the integration step.</p>
<p>In order to test for relative ordering, we need at least 2 markers in common (bridge markers) between each pair of homologue maps (Figure 1). Of course, more than 2 markers in common is better (and desirable for map integration) but we will consider it as the absolute minimum. A simple test to see whether the order is consistent as in Figure 1.a is then to perform a correlation on the marker positions, and see whether the result is positive or negative.</p>
<div class="figure">
<img src="" alt="Figure 1. a. Consistent marker ordering; b. inconsistent marker order" />
<p class="caption">Figure 1. a. Consistent marker ordering; b. inconsistent marker order</p>
</div>
<p>Both the correlation and integration steps have been put together into a single function, <code>orient_and_merge_maps</code>. <code>LPmerge</code> produces a number of integrated maps, based on an input parameter called the maximum interval size. A maximum interval size of one means that only adjacent markers on the underlying homologue maps are considered in the estimation of root mean square error (RMSE) in the final consensus map; a higher interval size means non-adjacent markers are also considered. The default setting is 4 – we recommend increasing this to 8. More details can be found in the literature on the subject <span class="citation">[3]</span>. The “best” map is generally the one with the minimum RMSE (but you may also like to select the integrated map with the smallest difference between the length of the consensus map and the mean map length of the homologue maps). Integrating 8 homologues of a single linkage group would look like so:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">integrated_map_LG2 &lt;-<span class="st"> </span><span class="kw">orient_and_merge_maps</span>(maplist_P1[[<span class="st">&quot;LG2&quot;</span>]], 
                                           maplist_P2[[<span class="st">&quot;LG2&quot;</span>]],
                                           <span class="dt">connection_threshold =</span> <span class="dv">3</span>,
                                           <span class="dt">plot_graph=</span><span class="ot">TRUE</span>)</code></pre></div>
<p>A code to provide a list of integrated maps of all linkage groups would look like this:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">integrated_maps &lt;-<span class="st"> </span><span class="kw">list</span>()
for(i in <span class="dv">1</span>:<span class="dv">5</span>){
  integrated_maps[[<span class="kw">paste0</span>(<span class="st">&quot;LG&quot;</span>, i)]] &lt;-<span class="st"> </span><span class="kw">orient_and_merge_maps</span>(maplist_P1[[<span class="kw">paste0</span>(<span class="st">&quot;LG&quot;</span>, i)]], 
                                                             maplist_P2[[<span class="kw">paste0</span>(<span class="st">&quot;LG&quot;</span>, i)]],
                                                             <span class="dt">LPmerge_interval =</span> <span class="dv">2</span>,
                                                             <span class="dt">connection_threshold =</span> <span class="dv">3</span>,
                                                             <span class="dt">plot_graph=</span><span class="ot">TRUE</span>)
}</code></pre></div>
<p><a name="sect15"></a></p>
</div>
<div id="plotting-a-map" class="section level2">
<h2>15. Plotting a map</h2>
<p>A well-know software program for plotting maps is MapChart <span class="citation">[4]</span>. <code>polymapR</code> can write out MapChart compatible files (.mct) using the function <code>write.mct</code>. There are many plotting options available in MapChart, here we only currently incorporate a single option, namely <code>showMarkerNames</code>, which by default is <code>FALSE</code> in anticipation of high-density maps. Further formatting can be achieved within the MapChart environment itself.</p>
<p>However, <code>polymapR</code> also comes with its own simple built-in function to plot maps, namely the <code>plot_map</code> function:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot_map</span>(<span class="dt">maplist =</span> integrated.maplist)</code></pre></div>
<p><img src="" /><!-- --></p>
<p>We might also want to visualise the integrated map output using the <code>plot_phased_maplist</code> function (but here we limit it to a single linkage group). Note that we need to have first used the <code>create_phased_maplist</code>, described earlier.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot_phased_maplist</span>(<span class="dt">phased.maplist =</span> phased.maplist[<span class="dv">1</span>], <span class="co">#Can plot full list also, remove &quot;[1]&quot;</span>
                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                    <span class="dt">cols =</span> <span class="kw">c</span>(<span class="st">&quot;black&quot;</span>,<span class="st">&quot;grey50&quot;</span>,<span class="st">&quot;grey50&quot;</span>))</code></pre></div>
<p><img src="" /><!-- --></p>
<p>This function can also visualise phased hexaploid maps, or triploid maps (tetraploid x diploid) if <code>ploidy2</code> is specified.</p>
<p><a name="sect16"></a></p>
</div>
<div id="evaluating-map-quality" class="section level2">
<h2>16. Evaluating map quality</h2>
<p><code>MDSMap</code> produces output that can be examined each time a map is produced, namely the principal coordinate plots with the principal curve (on which the map is based) - highlighting possible outlying markers that do not associate well with the rest of the markers, and a plot of the nearest-neighbour fit (nnfit) for the markers along each linkge group, giving an indication of markers with high stress. Both of these can help the user to evaluate the quality of the map (for more details, we recommend reading the <code>MDSMap</code> publication <span class="citation">[2]</span>) and often results in a number of rounds of mapping, where problematic markers are identified and removed followed by subsequent remapping and re-evaluation. Once the maps have been produced it is also an idea to get a general overview of the map quality using the <code>check_map</code> function. Note that by default this function expects lists as input arguments, enabling diagnostics to be run over all maps. Here we run it for a single linkage group, LG1:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">check_map</span>(<span class="dt">linkage_list =</span> <span class="kw">list</span>(linkages$LG1), <span class="dt">maplist =</span> <span class="kw">list</span>(integrated.maplist$LG1))</code></pre></div>
<p><img src="" width="600px" /></p>
<p>For each map produced, there are 3 diagnostic plots. The first of these shows the differences between the pairwise estimates of recombination frequency and the effective estimate of recombination frequency based on the map itself (the multi-point estimate of the recombination frequency) . These differences are compared to the LOD score for each pairwise estimate. An overall estimate of the weighted root mean square error of these differences (weighted by the LOD scores) is printed on the plot. In general, we expect that if the difference between the expected and realised recombination frequency is high, the LOD should be low. If this is not the case, something went wrong in the final mapping step. The second plot shows the comparison between a marker’s map position and the recombination frequency estimates to all other markers. In gerenal we expect that small recombination frequency estimates come from nearby markers and should therefore be concentrated around the diagonal, which are shown as light green regions. Areas of light green off the diagonal suggest a poor map order. The final plot is similar to the previous except that LOD values are shown in place of recombination frequencies. Again, we expect a concentration of higher LOD scores around the diagonal. By default, LOD scores less than 5 are not shown to make things clearer, but this value can be varied using the <code>lod.thresh</code> argument.</p>
<p><a name="sect17"></a></p>
</div>
<div id="preferential-pairing" class="section level2">
<h2>17. Preferential pairing</h2>
<p>Up to now, we have been working under the assumption of random bivalent pairing in the parents. However, it has been shown that in certain species the pairing behaviour is neither completely random nor completely preferential (as we have in allopolyploids) but something in-between, a condition termed segmental allopolyploidy <span class="citation">[5]</span>. A comprehensive study on this topic was recently performed in rose <span class="citation">[6]</span>, which found evidence of disomic behaviour in certain regions of the genome, with tetrasomic behaviour everywhere else. This sort of mixed inheritance can complicate the analysis and in cases where this effect is particularly pronounced, it is probably unwise to ignore <span class="citation">[7]</span>.</p>
<p>The <code>polymapR</code> package can currently accomodate preferential pairing in tetraploid species only, by correcting the recombination frequency estimates once the level of preferential pairing in a parental chromosome is known. In order to determine whether this is necessary, we can run some diagnostic tests using the marker data after an initial map has been produced which can inform our choice regarding whether to include a preferential pairing parameter or not (in a subsequent re-analysis).</p>
<p>We first test for preferential pairing using the function <code>test_prefpairing</code>, which re-examines closely-linked repulsion-phase simplex x nulliplex marker pairs. The method used here is described in Bourke et al.<span class="citation">[6]</span>, which is a development on ideas originally described by Qu and Hancock<span class="citation">[8]</span>. We look for signatures of preferential pairing in repulsion-phase pairs which map to the same locus. Defining a minimum distance by which to consider markers essentially at the same locus is somewhat arbitrary - in our study of rose, we used pairs of duplicated individuals to determine an approximate error rate which, along with considerations of the size of the population and rate of missing values, led us to conclude that distances less than ~1 cM were close enough to be informative. By default, we use an even more stringent value of 0.5 cM, which limits the number of marker pairs used in the comparison. This can be increased or decreased as required using the <code>min_cM</code> argument, as we show here:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">P1.prefPairing &lt;-<span class="st"> </span><span class="kw">test_prefpairing</span>(<span class="dt">dosage_matrix =</span> ALL_dosages, 
                                   <span class="dt">maplist =</span> integrated.maplist, 
                                   <span class="dt">LG_hom_stack =</span> LGHomDf_P1_1, 
                                   <span class="dt">min_cM =</span> <span class="dv">1</span>,
                                   <span class="dt">ploidy =</span> <span class="dv">4</span>)

<span class="kw">head</span>(P1.prefPairing)</code></pre></div>
<p>It is then a matter of checking the significance and deciding whether it is worthwhile to re-run the mapping. The column <code>P_value.adj</code> gives the adjusted P-value from a Binomial test to see whether the disomic repulsion estimate of recombination frequency differs significantly from <span class="math inline">\(\frac{1}{3}\)</span>. In cases where this is significant (with <code>P_value.adj</code> &lt; 0.01, say), we might have cause for concern. The output of <code>test_prefpairing</code> also lists the markers, their positions and on which homologues they reside. There is also an estimate of the pairwise preferential pairing parameter <code>pref.p</code>. If there is strong evidence of preferential pairing from multiple marker pairs, the approach we recommend is to take the average of the <code>pref.p</code> values for a particular chromosomal linkage group, which becomes one of the parameters <code>p1</code> or <code>p2</code> in the <code>linkage</code> function.</p>
<p>Here, we have purely random pairing in our parents, so the estimates of <code>pref.p</code> are nonsensical - in general, we have defined the preferential parameter <code>p</code> such that <span class="math inline">\(0 &lt; p &lt; \frac{2}{3}\)</span>. Supposing this were not the case, and we wanted to estimate the strength of preferential pairing on LG1 of parent 1, we would do something like the following:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">mean</span>(P1.prefPairing[P1.prefPairing$LG_a ==<span class="st"> </span><span class="dv">1</span>,]$pref.p)</code></pre></div>
<p>Finally, to re-run the linkage analysis with this estimate for parent 1 of linkage group 1, we would first have to subset our marker dosage data corresponding to that linkage group only (since we do not need to re-analyse any other chromosome!), and then re-run the linkage analysis. Note that in cases with extreme levels of preferential pairing, the clustering into separate homologues can also be affected. If we had estimated a value of 0.25 for the preferential pairing parameter (instead of -0.56, which is nonsense) we could proceed as follows:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">lg1_markers &lt;-<span class="st"> </span><span class="kw">unique</span>(<span class="kw">c</span>(<span class="kw">rownames</span>(marker_assignments_P1[marker_assignments_P1[,<span class="st">&quot;Assigned_LG&quot;</span>] ==<span class="st"> </span><span class="dv">1</span>,]),
                        <span class="kw">rownames</span>(marker_assignments_P2[marker_assignments_P2[,<span class="st">&quot;Assigned_LG&quot;</span>] ==<span class="st"> </span><span class="dv">1</span>,])))

filtered_data_lg1 &lt;-<span class="st"> </span>filtered_data[<span class="kw">rownames</span>(filtered_data) %in%<span class="st"> </span>lg1_markers,]

all_linkages_list_P1_lg1 &lt;-<span class="st"> </span><span class="kw">finish_linkage_analysis</span>(<span class="dt">marker_assignment =</span> marker_assignments$P1[<span class="kw">rownames</span>(marker_assignments$P1) %in%<span class="st"> </span>lg1_markers,],
                                                    <span class="dt">dosage_matrix =</span> filtered_data_lg1,
                                                    <span class="dt">target_parent =</span> <span class="st">&quot;P1&quot;</span>,
                                                    <span class="dt">other_parent =</span> <span class="st">&quot;P2&quot;</span>,
                                                    <span class="dt">convert_palindrome_markers =</span> <span class="ot">FALSE</span>,
                                                    <span class="dt">ploidy =</span> <span class="dv">4</span>,
                                                    <span class="dt">pairing =</span> <span class="st">&quot;preferential&quot;</span>,
                                                    <span class="dt">prefPars =</span> <span class="kw">c</span>(<span class="fl">0.25</span>,<span class="dv">0</span>),
                                                    <span class="dt">LG_number =</span> <span class="dv">1</span> <span class="co">#interested in just 1 chm.</span>
)</code></pre></div>
<p><a name="sect18"></a></p>
</div>
<div id="concluding-remarks" class="section level2">
<h2>18. Concluding remarks</h2>
<p>We hope you have been able to successfully create an integrated map using the methods described here. We anticipate <code>polymapR</code> will improve over the coming years (in terms of applicability as well as finding and removing bugs). Feedback on performance including bug reporting is most welcome; please refer to the package information on CRAN for details on who to contact regarding maintenance issues.</p>
<p><a name="sect19"></a></p>
</div>
<div id="references" class="section level2 unnumbered">
<h2>19. References</h2>
<div id="refs" class="references">
<div id="ref-Voorrips2011">
<p>1. Voorrips RE, Gort G, Vosman B: <strong>Genotype calling in tetraploid species from bi-allelic marker data using mixture models</strong>. <em>BMC bioinformatics</em> 2011, <strong>12</strong>:172.</p>
</div>
<div id="ref-Hackett2016">
<p>2. Preedy KF, Hackett CA: <strong>A rapid marker ordering approach for high-density genetic linkage maps in experimental autotetraploid populations using multidimensional scaling</strong>. <em>Theoretical and Applied Genetics</em> 2016, <strong>129</strong>:2117–2132.</p>
</div>
<div id="ref-Endelman2014">
<p>3. Endelman JB, Plomion C: <strong>LPmerge: An R package for merging genetic maps by linear programming</strong>. <em>Bioinformatics</em> 2014, <strong>30</strong>:1623–1624.</p>
</div>
<div id="ref-Voorrips2002">
<p>4. Voorrips R: <strong>MapChart: software for the graphical presentation of linkage maps and QTLs</strong>. <em>Journal of heredity</em> 2002, <strong>93</strong>:77–78.</p>
</div>
<div id="ref-Stebbins1947">
<p>5. Stebbins GL: <strong>Types of polyploids: their classification and significance</strong>. <em>Advances in genetics</em> 1947, <strong>1</strong>:403–429.</p>
</div>
<div id="ref-Bourke2017">
<p>6. Bourke PM et al: <strong>Partial preferential chromosome pairing is genotype dependent in tetraploid rose</strong>. <em>The Plant Journal</em> 2017, <strong>90</strong>:330–343.</p>
</div>
<div id="ref-Bourke2016">
<p>7. Bourke PM, Voorrips RE, Kranenburg T, Jansen J, Visser RGF, Maliepaard C: <strong>Integrating haplotype-specific linkage maps in tetraploid species using SNP markers</strong>. <em>Theoretical and Applied Genetics</em> 2016, <strong>129</strong>:2211–2226.</p>
</div>
<div id="ref-Qu2001">
<p>8. Qu L, Hancock J: <strong>Detecting and mapping repulsion-phase linkage in polyploids with polysomic inheritance</strong>. <em>Theoretical and Applied Genetics</em> 2001, <strong>103</strong>:136–143.</p>
</div>
</div>
</div>



<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>
